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ABSTRACT 


The  use  of  cellular  phones  has  become  widespread.  It  is  predicted  that  cellular 
phone  use  will  soon  become  almost  as  high  as  normal  land-line  telephones.  The  ability  to 
track  these  phones  has  obvious  advantages  for  intelligence  gathering.  The  problem  with 
tracking  these  phones  is  that  their  use  is  intermittent  and  they  are  very  low  power 
emitters.  The  use  of  Unmaimed  Aerial  Vehicles  (UAV’s)  could  help  in  the  detection  of 
these  signals.  They  would  extend  the  capabilities  of  a  ground  unit.  Once  these  signals 
have  been  detected,  the  best  algorithm  for  tracking  them  is  the  Kalman  filter  because  of 
the  intermittent  and  noisy  nature  of  the  received  signals.  The  extended  Kalman  filter  is 
used  because  of  the  nonlinearities  present  in  the  system.  The  exploitation  of  cellular 
signals  by  using  UAV's  and  the  extended  Kalman  filter  is  an  important  firework  for 
future  use  of  UAV's  in  unconventional  ways.  This  thesis  explores  the  use  of  UAV’s  to 
exploit  cellular  emissions  and  locate  the  caller.  The  time  difference  of  arrival  (TDOA)  of 
emissions  is  the  main  method  of  tracking  using  the  Kalman  filter.  Further  research  and 
directions  of  interest  will  be  proposed. 
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I.  INTRODUCTION 


A.  BACKGROUND 

Cellular  phones  are  becoming  an  everyday  appliance  to  many  people.  The 
widespread  use  of  these  phones  has  sparked  interest  in  the  intelligence  community, 
because  cellular  phone  emission  characteristics  are  such  that  the  signals  can  be  tracked. 
By  locating  a  certain  phone,  the  operator's  position  is  likewise  known.  Being  able  to 
locate  individuals  to  a  certain  degree  of  accuracy  could  be  very  beneficial.  Currently 
there  are  few  ways  to  track  these  phones  in  a  passive  way.  This  work  will  focus  upon  the 
problem  of  finding  the  position  of  a  cellular  phone. 

The  UAV  (Unmanned  Aerial  Vehicle)  has  been  very  useful  in  providing 
battlefield  surveillance,  target  information,  and  naval  gunfire  spotting.  The  use  of  UAV’s 
allows  the  intelligence  field  the  flexibility  to  obtain  information  that  is  too  risky  to  obtain 
with  manned  aircraft  or  impossible  to  get  with  a  satellite.  The  UAV  provides  some 
measure  of  stealth  since  its  small  size  makes  it  a  hard  target  for  radar  to  detect.  The  small 
size  also  is  of  benefit  when  the  UAV  becomes  a  target,  making  it  harder  to  hit. 

The  UAV  could  be  a  valuable  asset  in  tracking  cellular  phones.  A  small  UAV 
could  be  carried  in  the  back  of  a  truck  or  in  a  small  boat.  The  UAV,  when  fitted  with  a 
simple  omni-directional  antenna,  would  be  able  to' receive  cellular  transmissions.  By 
fitting  the  ground/water  vehicle  with  a  similar  antenna,  one  could  measure  the  time 
difference  of  arrival  (TDOA)  for  the  reception  of  a  cellular  emission  between  the  base 
receiver  and  the  UAV.  This  TDOA  is  used  to  estimate  the  position  of  the  phone.  The 
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addition  of  the  UAV  to  the  intelligence  collection  vehicle  gives  the  operator  a  second 
receiver,  which  makes  the  use  of  these  calculations  possible. 

B.  CELLULAR  PHONES 

The  excellent  work  of  Lee  is  used  as  the  main  reference  when  dealing  with 
cellular  phones  [Ref.  1].  Some  consider  it  the  bible  of  cellular  phones. 

In  general,  cellular  systems  are  classified  by  operating  frequencies  and  channel 
spacing.  All  systems  operate  at  either  450  MHz  or  800  MHz  and  have  channel 
bandwidths  of  30, 25,  or  20  kHz.  There  are  also  two  basic  types  of  transmissions  of  voice 
data  for  these  systems,  analog  and  digital.  Today,  digital  systems  are  just  being 
introduced  in  the  commercial  market.  Digital  systems  allow  a  greater  number  of  callers 
and  have  lower  interference  problems  than  the  analog  systems. 

The  system  of  choice  in  North  America  is  the  Advanced  Mobile  Phone  Service 
(AMPS),  which  was  first  placed  into  service  in  1978.  The  rest  of  the  world  uses  various 
systems  that  are  similar  to  AMPS,  but  not  compatible.  This  work  will  focus  on  gathering 
intelligence  from  the  AMPS  cellular  system.  The  methods  employed  can  be  modified  to 
work  with  the  other  cellular  systems  of  the  world. 

C.  CELLULAR  PHONE  EMISSIONS 

Cellular  phone  emissions  are  unique  in  many  ways.  The  cellular  phone  system  is 
based  upon  the  phone  cell  and  its  associated  firequency  reuse  pattern.  A  cell  is  a 
geographic  area  which  is  assigned  certain  frequencies  of  the  total  spectrum.  The  size  of 
the  area  depends  upon  the  expected  traffic  in  the  cell.  An  area  of  operation  is  broken 
down  into  a  number  of  these  cells  with  each  cell  covering  a  unique  portion  of  the 
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operating  area,  as  illustrated  in  Figure  1.  This  is  an  ideal  representation,  as  most  cells 
will  not  have  a  circular  shape  and  will  be  of  different  sizes. 


Figure  1.  Cell  Structure. 

The  cellular  structure  shown  is  based  upon  a  seven  cell  reuse  pattern.  What  this 
means  is  that  each  cell  of  a  reuse  pattern  is  assigned  a  certain  segment  of  the  total  cellular 
channels  available.  A  seven  cell  reuse  pattern  means  that  the  total  number  of  voice 
channels  is  broken  up  seven  ways  and  split  among  the  cells.  Cells  with  the  same  channel 
number  use  the  same  voice  channels.  It  is  important  to  separate  cells  using  the  same 
frequencies  by  a  distance  that  will  insure  that  transmissions  will  not  interfere  with  each 
other.  This  interference  is  called  cross-talk.  To  prevent  cross-talk,  the  cells  are  kept 
small  and  the  power  outputs  are  kept  low.  The  AMPS  system  specifies  that  the  cell 
radius  be  from  2  km  to  20  km.  The  use  of  smaller  cells  allows  a  greater  number  of  callers 
to  be  serviced  in  the  same  size  subscription  area.  This  cell  structure  causes  a  problem 
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when  trying  to  track  one  phone,  since  there  may  be  a  number  of  phones  using  the  same 
frequencies  in  a  given  area,  and  the  phone  will  switch  frequencies  as  it  travels  across  cell 
boundaries.  This  switching  of  frequencies  as  the  phone  moves  is  called  hand-off. 

The  AMPS  system  is  an  analog  system.  This  means  that  voice  information  is 
transmined  using  frequency  modulation.  Each  channel  is  split  into  a  forward  voice 
channel  (FVC)  and  a  reverse  voice  channel  (RVC).  The  forward  voice  channel  transmits 
from  the  base  station  to  the  mobile  station  and  the  reverse  voice  channel  from  the  mobile 
station  to  the  base  station.  The  forward  and  reverse  channels  have  a  bandwidth  of  30  kHz 
and  are  separated  from  each  other  by  45  kHz.  In  the  AMPS  system,  forward 
transmissions  (station  to  phone)  are  from  870  to  890  kHz  and  the  reverse  transmissions 
(phone  to  station)  are  from  825  to  845  kHz  as  can  be  seen  in  Figure  2.  The  spectrum  is 
up  of  a  total  of  666  chaimels,  t^ch  are  divided  into  two  blocks,  A  and  B,  so  that 
fftmpfting  cellular  companies  can  each  use  half  of  the  total  spectrum  as  delineated  by  the 
FCC.  The  chaimels  313  to  354  are  reserved  as  control  channels,  and  the  rest  are  used  for 
voice  transmission,  giving  each  carrier  a  total  of  21  setup  control  chaimels  and  312  voice 
oti«rtwii>u  On  26  July  1986,  the  FCC  increased  available  spectrum  by  10  MHz,  which 
added  166  more  voice  channels.  The  21  control  channels  are  frequency  shift  keyed 
(FSK)  digital  signals. 
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Figure  2.  Cellular  Channel  Spectrum  [After  Ref.  1]. 


The  control  channels  perform  a  number  of  tasks,  including  channel  setup.  They 
can  be  access  channels,  used  to  setup  mobile  originated  calls,  or  paging  channels,  used 
for  land  originated  calls.  These  signals  are  digital,  and  are  modulated  with  frequency 
shift  keying  (FSK).  Digital  messages  are  sent  in  a  Burst  mode  with  a  10  kbps 
transmission  rate.  There  are  also  two  analog  control  signals,  the  supervisory  audio  tone 
(SAT)  and  the  signaling  tone  (ST). 

The  SAT  is  of  specific  interest  in  this  thesis.  It  is  a  tone  that  is  added  to  each 
FVC  when  an  active  link  is  made,  and  is  placed  outside  of  the  bandwidth  that  transmits 
the  voice  information,  but  within  the  channel  bandwidth.  The  mobile  phone  detects  this 
SAT,  filters  it,  and  then  returns  it  on  the  RVC.  This  was  originally  designed  to  allow  the 
base  station  to  determine  the  phase  difference  when  the  SAT  returned.  This  phase 
difference  could  be  used  to  compute  the  total  round  trip  time  from  base  to  mobile  and 
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back  to  the*  base,  which  could  then  be  used  to  calculate  distance  from  the  base  station  to 


the  mobile  phone.  It  is  not  used  for  distance  calculations  now,  but  serves  other  control 
purposes.  Today  most  systems  determine  when  to  hand-off  to  the  next  cell  based  upon 
the  signal  strength  of  the  RVC  or  the  power  difference  of  the  mobile  signal  received  by 
adjacent  cell  sites. 

D.  CELLULAR  PHONE  EQUIPMENT 

A  typical  cellular  phone  area  is  divided  into  a  large  number  of  cells  which  use  a 
certain  frequency  reuse  scheme.  These  cells  may  not  be  the  same  size  or  shape, 
depending  upon  the  propagation  characteristics  of  the  emissions  &om  the  ceil  site 
antennae  and  the  expected  traffic.  Each  cell  site  consists  of  antennae,  radios,  control 
units,  and  data  terminals.  These  cell  sites  are  all  linked  in  various  ways  to  the  Mobile 
Telephone  Switching  Office  (MTSO).  It  is  here  that  all  the  switching  is  accomplished 
between  the  cellular  network  and  the  land  lines;  the  calls  are  routed  and  connected  to  the 
nationwide  telephone  network.  The  MTSO  is  the  heart  of  the  cellular  system  and 
provides  central  coordination  and  administration  for  the  system. 

Mobile  phones  have  become  very  small  in  the  last  couple  of  years.  The 
transmitting  power  is  relatively  low,  with  the  maximum  usually  4  watts.  The  phones  do 
not  always  transmit  at  the  maximum  power  level,  and  usually  transmit  just  the  minimum 
power  required  for  the  base  station  to  receive  the  signal.  This  is  done  to  reduce  the 
interference  between  phones.  In  the  future  digital  phones  will  become  even  more 
widespread,  as  they  are  better  suited  for  mobile  data  communications. 
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E.  PHONE  LOCATION  PROBLEM  DESCRIPTIONS 

Two  basic  approaches  are  used  to  estimate  the  location  of  a  cellular  phone  from 
TDOA  information.  In  the  first  approach  the  SAT  is  used  to  estimate  the  total  distance 
the  signal  has  traveled  from  the  base  station  to  the  phone  and  then  to  the  receiver.  An 
ellipse  of  estimated  location  is  then  calculated.  The  second  approach  does  not  involve  the 
emissions  from  the  base  station.  In  this  method  the  digital  Burst  emissions  from  the 
phone  are  marked  in  time.  The  TDOA  between  two  different  receivers  is  then  used  to 
estimate  position.  The  two  approaches  may  be  combined  if  both  forms  of  information  are 
available. 

For  both  of  these  methods  it  will  be  assumed  that  one  receiver  is  based  out  of  a 
stationary  van  (or  some  other  collection  vehicle).  In  addition  to  the  stationary  vehicle,  the 
user  will  have  at  his  disposal  one  UAV  that  will  be  able  to  receive  cellular  emissions.  The 
UAV  will  be  flown  in  a  certain  pattern  that  optimizes  the  estimation  algorithm. 

F.  SAT  TDOA  PROBLEM  DESCRIPTION 

The  first  TDOA  filtering  problem  assumes  that  a  call  is  in  progress  between  the 
desired  phone  and  a  cell  base  station.  The  detection  of  this  condition  involves  listening  in 
the  area  of  estimated  operation  and  waiting  for  the  control  signals  that  setup  the  call  and 
assign  the  channel  frequencies.  This  is  easy  if  the  call  is  originated  from  land  lines,  since 
the  cellular  system  "pages"  the  cellular  phone  by  broadcasting  a  setup  message  from 
every  cell  in  the  subscription  area.  When  the  phone  responds,  the  system  knows  which 
cell  to  use.  Calls  originating  from  a  cellular  phone  are  harder  to  detect.  The  unique 
identification  number  for  each  phone  will  be  the  flag  that  announces  that  a  channel 
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assignment  is  about  to  be  made.  Once  the  channel  is  known,  the  SAT  from  the  FVC  and 
RVC  can  be  easily  timed  in.  In  this  problem  the  following  assumptions  are  made: 

1 .  Enough  information  is  obtained  from  the  received  signals  to  filter  out 
other  emitters.  This  reduces  the  multi-emitter  problem  into  a  single 
emitter  problem. 

2.  The  receivers  are  referenced  to  a  conunon  time  base.  This  could  be 
done  by  using  the  clock  in  the  Global  Positioning  System  (GPS) 
system. 

3.  GPS  will  be  used  in  both  the  land  receiver  and  the  UAV.  The 
positions  of  the  receivers  will  be  assumed  to  be  without  error  since  the 
accuracy  of  the  GPS  system  is  better  than  our  desired  system  accuracy. 

4.  The  receivers  each  use  one  omnidirectional  antenna.  This  will  keep 
system  hardware  as  simple  as  possible. 

5.  The  data  link  between  the  UAV  and  the  intelligence  van  will  be 
assumed  to  be  fast  enough. 

6.  The  UAV  will  fly  some  set  pattern.  The  efficiencies  of  these  patterns 
will  be  discussed.  A  standard  pattern  will  allow  automatic  flight 
control  and  relieve  the  operator  from  an  additional  task. 

7.  The  phone  in  question  will  be  stationary  for  the  duration  of 
measurements.  This  will  eliminate  the  need  to  switch  frequencies  as 
hand-offs  occur. 

8.  The  exact  location  of  the  cellular  base  station  in  use  will  be  known 
accurately. 

The  ranges  in  question  here  are  limited  to  the  size  of  one  cell.  It  will  be  assumed 
that  the  cell  in  use  has  a  diameter  of  20  km. 

G.  BURST  TDOA  PROBLEM  DESCRIPTION 

The  Burst  TDOA  problem  estimates  the  location  of  the  emitter  by  measuring  the 
TDOA  of  an  emission  from  the  cellular  phone  between  the  intelligence  van  and  the  UAV. 
The  position  of  the  cell  base  station  is  not  needed  and  emissions  from  it  are  not  used.  The 
following  assumptions  are  made  in  developing  the  problem: 
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1.  Some  sort  of  DSP  method  will  be  used  to  detect  the  digital  FSK 
control  bursts  that  the  phone  periodically  emits. 

2.  The  receivers  are  referenced  to  a  common  time  base.  This  could  be 
done  by  using  the  clock  in  the  Global  Positioning  System  (GPS) 
system. 

3.  GPS  will  be  used  in  both  the  land  receiver  and  the  UAV.  The 
positions  of  the  receivers  will  be  assumed  to  be  without  error,  since  the 
accuracy  of  the  GPS  system  is  better  than  the  desired  system  accuracy. 

4.  The  receivers  each  use  one  omnidirectional  antenna.  This  will  keep 
system  hardware  as  simple  as  possible. 

5.  The  data  link  between  the  UAV  and  the  intelligence  van  will  be 
assumed  to  be  fast  enough  for  our  purposes. 

6.  The  phone  in  question  will  be  stationary  for  the  duration  of 
measurements.  This  will  eliminate  the  need  to  switch  frequencies  as 
hand-offs  occur. 

7.  The  UAV  will  fly  some  set  pattern.  The  efficiencies  of  these  patterns 
will  be  discussed.  A  standard  pattern  will  allow  automatic  flight 
control  and  relieve  the  operator  from  an  additional  task. 

8.  the  exact  location  of  the  cellular  base  station  in  use  will  not  be  known. 

This  algorithm  will  also  be  tested  for  the  area  of  one  cell,  a  20  km  region.  It  could 
be  used  for  a  much  larger  region  when  data  bursts  can  be  received  regardless  of  which 
cell  the  phone  is  in.  In  theory  this  is  possible,  but  it  must  be  remembered  that  there  will 
be  a  great  number  of  phones  using  the  same  control  channels.  A  large  portion  of  work, 
which  will  not  be  done  here,  will  be  required  to  be  able  to  find  the  desired  ID  number. 
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II.  TDOA  ESTIMATION 


A.  TIME  OF  ARRIVAL  ESTIMATION  AND  ERROR  MODELS 

The  Kalman  filter  is  used  to  estimate  the  position  of  the  cellular  phone  from  the 
TDOA  observations.  The  Kalman  filter  is  used  to  estimate  the  instantaneous  state  of  a 
system  perturbed  by  Gaussian  white  noise  by  using  measmements  that  are  related  to  the 
states,  but  are  corrupted  by  Gaussian  white  noise.  The  filter  is  a  statistically  optimal 
estimator.  In  this  case,  the  noise  of  interest  will  be  the  measurement  noise,  rather  than  the 
process  noise.  To  use  the  Kalman  filter  this  noise  will  be  assiuned  to  be  Gaussian  white. 

This  assumption  can  be  made  because  most  real  noise  is  in  fact  close  enough  to 
this  ideal  so  that  the  filter  algorithms  will  work.  If  the  actual  measurement  noise  is  not 
white,  but  is  instead  colored,  it  can  be  whitened  by  using  shaping  filters. 

B.  TIME  OF  ARRIVAL  MEASUREMENT 

TDOA  measurements  between  the  van  and  the  UAV  will  be  used  to  estimate  the 
location  of  the  phone.  Since  there  is  no  instantaneous  physical  link  between  the  two 
receivers,  it  is  important  that  the  clocks  in  each  are  synchronized  precisely.  To  minimize 
the  drift  that  will  occur  between  any  two  clocks,  they  are  synchronized  to  the  GPS 
internal  clock .  The  GPS  uses  this  clock  in  its  own  position  calculations.  Further  checks 
on  clock  synchronization  could  be  made  over  the  data  link  by  using  any  of  the  time  check 
algorithms  currently  in  use  in  computer  networks.  This  assiunes  that  the  link  between  the 
van  and  UAV  is  in  the  form  of  digital  communication.  This  link  could  be  formed  over 
some  sort  of  wireless  link,  fiber  optic  line,  or  even  stored  for  download  upon  landing  the 
UAV. 
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For  the  SAT  TDOA  observations  the  variations  encountered  in  the  GPS  clock 
synchronization  should  be  too  small  to  affect  the  system.  The  SAT  tone  is  typically  at  6 
kHz.  This  gives  a  peak  to  peak  time  of  166.6  microseconds.  The  GPS  clock 
synchronization  error  will  be  a  couple  of  orders  of  magnitude  below  this. 

The  Burst  TDOA  observation  accuracy  depends  upon  the  type  of  detection 
algorithm  used  to  mark  the  time  of  arrival  (TOA).  The  interference  may  be  so  great  and 
the  signal  so  short  that  either  receiver  or  both  may  miss  the  signal  altogether.  This  also 
depends  upon  what  part  of  the  cellular  emission  is  exploited.  This  algorithm  could 
trigger  the  emission  of  the  ID  code  from  the  cellular  phone,  which  would  relieve  the  great 
detection  burden,  but  assumes  that  the  phone  is  engaged  in  a  call.  In  the  future,  cellular 
phones  will  periodically  send  a  registration  ID  to  the  cell  site  so  that  the  MTSO  will  be 
able  to  track  phone  locations.  This  type  of  autonomous  registration  will  be  assumed  to  be 
what  the  Burst  TDOA  model  is  receiving. 

C.  SAT  TIME  OF  ARRIVAL  ESTIMATION 

The  SAT  is  a  6  kHz  tone  added  to  the  RVC  and  FVC.  Since  it  is  of  relatively  low 
frequency,  it  is  possible  to  mark  the  reception  time  from  the  peaks  of  the  signal.  As  long 
as  the  actual  time  difference  is  not  larger  than  the  peak  to  peak  time,  there  will  be  no 
confusion  as  to  which  received  peaks  match  the  two  receivers.  The  period  of  the  signal  is 
1 66.7  microseconds.  Multiplying  this  by  the  speed  of  light  yields: 

3  X  10*  X  166.67  X  10-^  =  50  X 10^  (1) 
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Thus,  as  long  as  the  difference  between  the  two  receivers  is  less  than  50  km,  this 
signal  is  used.  This  does  not  mean  that  the  range  of  detection  is  limited  to  50  km,  only 
that  the  separation  between  the  receivers  should  be  less  than  50  km. 

D.  BURST  TIME  OF  ARRIVAL  ESTIMATION 

The  Burst  TDOA  model  receives  its  time  of  arrival  information  after  the  signal 
processing  algorithm  determines  that  it  has  in  fact  received  the  desired  ID.  Since  the 
signal  is  a  digital  signal  that  manifests  itself  as  a  FSK  modulated  signal,  any  one  of  the 
established  computer  communication  algorithms  could  be  modified  for  use.  Of 
importance  is  the  point  where  the  algorithm  triggers  the  time  measurement.  It  would  be 
preferable  to  have  it  trigger  on  the  first  received  bit.  This  would  require  some  type  of 
memory  that  might  be  prohibitive.  Therefore,  the  triggering  will  occur  on  the  last  bit. 

As  stated  earlier,  since  the  signal  to  noise  ratio  will  be  low  and  the  algorithm  may 
not  recognize  the  signal,  the  reception  of  TOA  measurements  will  not  occur  at  a  uniform 
rate,  but  will  be  modeled  as  a  random  process. 

E.  SAT  AND  BURST  TDOA 

The  TOA  for  these  signals  is  approximated  as  a  Gaussian  random  variable.  It  is 
assumed  that  each  TOA  is  independent.  The  TDOA  is  calculated  by  using  the  following 
equation: 

TDOA=z^-z,^,  (2) 

where  =  TOA  of  the  emission  at  the  van,  and  -  TOA  of  the  emission  at  the 
UAV. 
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The  TOA  can  be  expressed  as  the  sum  of  a  true  TOA  and  a  Gaussian  white  noise 
random  variable  such  as 

z=Zo  +  v,  (3) 

where  z  =  the  observed  TOA,  Zq  =  the  true  TOA,  and  v  =  the  noise  present  in  the 
measurement. 

The  mean  of  the  noise  is  zero  and  has  a  variance  of  V.  The  expected  value  of 
E[VvanVuav]  IS  zero,  SO  the  noise  components  are  uncorrelated.  The  mean  of  the  TDOA 
estimates  is 

E[zv«  -  Zuav]  =  E[Zo,„  -  Zo,„v]+  E[Vv«]-  E[v„,v]  =  Zova  -  Zou»  •  W 

It  can  be  seen  that  the  expected  value  of  TDOA  is  the  difference  between  the  true  TOA. 
The  variance  will  be 

<jU.=E^TDOA-E[z,.-z„,])‘]=Efv_-v^y].  (5) 

The  two  noise  sequences  are  independent, 

E[v^v^J=0  (6) 

and 

=  E[vL  +  +  (7) 

Thus,  the  TDOA  is  a  random  variable  with  a  mean  value  equal  to  the  true  TDOA  and  a 
variance  that  is  equal  to  the  sum  of  the  variances  of  the  TOA  measurements. 
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III.  PHONE  LOCATION 


The  TDOA  is  the  time  difference  that  a  signal  or  series  of  signals  have  when  they 
are  received  by  receivers  that  are  separated  by  some  distance.  By  analyzing  the  geometry 
of  the  problems  the  location  of  the  phone  is  estimated.  This  location  estimate  will  not  be 
a  precise  point  when  only  one  receiver  is  used,  but  will  be  an  ellipse  or  locus  of  the 
possible  positions.  The  use  of  two  receivers  allows  the  pinpointing  of  the  phone  position 
by  looking  at  the  intersection  of  the  two  possible  position  loci  (ellipses). 

A.  USING  SAT  INFORMATION 
1.  Problem  Geometry 

In  this  problem,  the  position  of  the  cell  transmitter  is  assumed  to  be  known.  The 
receiver  measures  the  TDOA  between  the  SAT  from  the  base  station  on  the  FVC  and  the 
SAT  from  the  mobile  phone  on  the  RVC.  The  geometry  for  a  single  receiver  can  he  seen 
in  Figure  3. 


Figure  3.  Geometry  of  SAT  Problem. 
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The  originating  SAT  travels  from  the  base  station  (B)  to  the  phone  (P).  The 
phone  processes  the  signal  and  sends  it  back.  The  receiver  measures  the  TDOA  between 
the  emission  from  the  base  station  and  the  reply  SAT  from  the  phone.  In  this  problem  the 
following  values  are  known: 

d  =  the  distance  between  the  base  station  and  the  receiver, 
xb,  yb  =  the  position  of  the  base  station, 
xr,  yr  =  the  position  of  the  receiver, 

Tb  =  the  TOA  for  the  emission  from  the  base  station, 

Tp  =  the  TOA  for  the  emission  from  the  phone,  and 
c  =  the  speed  of  transmission,  taken  to  be  the  speed  of  light. 

The  following  values  are  unknown: 

xp,  yp  =  the  position  of  the  phone, 

Rbp  =  the  distance  between  the  base  station  and  the  phone, 

RpR  =  the  distance  between  the  phone  and  the  receiver. 

To  =  the  time  the  base  station  originates  the  transmission. 

The  time  of  origination  can  be  calculated  as 

To=Tb-%.  (8) 

The  time  it  takes  to  travel  from  the  base  station  to  the  phone  to  the  receiver,  not  counting 
processing  time  at  the  phone,  is 


Substituting  (8)  into  (9)  and  solving  for  the  distance 

R.,  +  Rp.  =  (t,  -(t,  -|]]c  =  [t,  -  T, +|)c,  (10) 

which  further  reduces  to 

Rbp  +  Rpr  =  TDOAxc  +  d.  (11) 

Equation  (11)  describes  an  ellipse  as  shown  in  Figure  3.  A  formula  can  be 
constructed  to  determine  the  coordinates  of  the  phone  based  upon  the  base  station 
position,  the  receiver  position,  and  the  TDOA. 

Of  more  importance  in  implementing  the  Kalman  filter  is  the  estimation  of  the 
TDOA  based  upon  the  known  position  of  the  base  station  and  receiver  and  the  estimated 
position  of  the  phone.  Solving  (11)  for  TDOA  yields: 

TDOA = +  >/(^p "  + (^  -  ^  >  (12) 

where 


2.  Loci  of  Constant  TDOA 

The  locus  of  constant  TDOA  is  an  ellipse  that  has  its  foci  at  the  base  station  and 
the  receiver.  The  phone  can  be  in  an  infinite  number  of  possible  locations  on  the  ellipse 
that  results  in  the  specific  TDOA.  The  problem  is  to  now  find  the  specific  location  on  the 
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locus.  One  solution  is  to  give  the  receiver  a  directional  antenna  so  that  a  bearing  can  be 
measured.  Using  the  ellipse  locus  and  the  bearing,  the  position  can  be  estimated. 

As  noted  before,  this  problem  utilizes  two  receivers  that  do  not  have  bearing 
capability.  By  looking  at  where  the  loci  of  each  receiver  intersect,  an  estimate  of  position 
can  be  made.  Notice  that  having  two  receivers  does  not  entirely  remove  the  ambiguities, 
since  they  intersect  in  two  places  (see  Figure  4).  To  pinpoint  which  of  these  intersections 
is  valid,  a  third  receiver  could  be  used.  Since  the  UAV  is  mobile,  each  measurement 
taken  as  it  moves  is  considered  to  be  a  different  measurement  from  a  separate  receiver. 
The  two  receivers  then  look  like  countless  receivers. 


UAV  loci 


Figure  4.  Loci  of  Constant  TDOA  Intersection. 
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B.  ORTHOGONALITY  OF  TDOA  OBSERVATIONS 

The  Kalman  filter  estimation  depends  a  great  deal  upon  how  orthogonal  the  two 
loci  of  positions  are.  The  more  orthogonal  the  intersection,  the  less  error  inherent  in  the 
system.  On  the  other  hand,  the  more  parallel  the  loci  are,  the  worse  the  estimation  will 
be.  This  inherent  error  can  be  shown  by  calculating  error  ellipsoids.  This  is  done  later  in 
this  work. 

It  is ‘desirable  to  have  the  loci  orthogonal  to  each  other.  Since  one  of  the  receivers 
(the  UAV)  is  mobile,  this  orthogonality  is  easily  obtained  through  maneuvers. 

C.  USING  ENCODED  BURST  INFORMATION 

1.  Problem  Geometry 

This  problem  does  not -consider  the  location  of  the  base  station,  since  only 
transmissions  from  the  phone  will  be  processed.  The  TDOA  will  consist  of  the  time 
difference  between  the  reception  times  of  a  phone’s  emissions  between  two  receivers. 
The  TDOA  is  a  fimction  of  the  distance  between  the  receivers.  The  problem  is  illustrated 
in  Figure  5. 
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Figures.  Pulse TDOA Geometry. 
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Xvan.yvan  =  the  position  of  the  van, 

Xuav.yuav  =  the  position  of  the  UAV,  and 
c  =  the  speed  of  light. 

2.  Loci  of  Constant  TDOA 

The  loci  of  constant  TDOA  are  hyperbolas  with  each  receiver  acting  as  a  focus 
(see  Figure  6).  There  is  only  one  locus  that  passes  through  the  middle  of  the  two 
receivers.  Here  again,  there  are  an  infinite  number  of  positions  where  the  phone  could  be. 
In  order  to  resolve  this  conflict,  three  receivers  would  allow  the  phone  to  be  pinpointed 
by  finding  where  the  three  loci  of  position,  which  pass  between  each  pair,  intersect.  Only 
two  receivers  are  used  in  this  problem,  but  the  mobility  of  the  UAV  allows  the  simulation 
of  more  than  two  receivers  by  taking  readings  at  different  positions  for  the  UAV. 


Loci  of  constant  TDOA 


Figure  6.  Loci  for  Burst  Data. 
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3.  Orthogonality  of  the  TDOA  Observations 

Kalman  filter  accuracy  depends  on  how  orthogonal  the  loci  of  position  are  at  the 
point  of  intersection.  As  mentioned  earlier,  the  more  orthogonal,  the  lower  the  expected  ' 
error.  This  will  be  shown  using  error  ellipses. 

4.  Moving  Phone 

In  this  work,  the  phone  is  assumed  to  be  stationary.  If  the  phone  moves,  a 
dynamic  state  is  assumed.  The  Kalman  filter  may  still  be  able  to  track  the  phone. 
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IV.  KALMAN  FILTERING 


A.  THE  EXTENDED  KALMAN  FILTER 

As  will  be  demonstrated  later,  the  plant  equations  employed  are  nonlinear.  For 
this  reason  the  normal  Kalman  filter  cannot  be  used,  and  the  extended  Kalman  filter  must 
be  substituted.  This  involves  the  use  of  partial  derivatives  as  linear  approximations  of  the 
nonlinear  models.  These  partial  derivatives  are  evaluated  at  the  estimated  values  of  the 
state  variables  and  are  used  to  calculate  the  Kalman  filter  gains. 

For  the  model  the  equations  of  state  are 

x(k  +  l)  =  f(x(k),w(k),k),  (17) 

where:  x(k)  =  the  state  of  the  system, 
w(k)  =  the  plant  driving  noise, 
f(k)  =  a  function  of  the  states,  noise,  or  k. 

The  observation  equations  are 

z(k)  =  h(x(k),v(k),k),  (18) 

where:  z(k)  =  the  system  measurements, 
v(k)  =  the  measurement  noise, 
h(k)  =  a  function  of  the  states,  noise,  or  k. 

Notice  that  both  f(k)  and  h(k)  can  be  nonlinear  functions.  In  both  of  the  problems 
h(k)  is  nonlinear.  The  two  noises,  w(k)  and  v(k)  are  assumed  to  be  white  noises  with  zero 
means  and  covariance's  of  W  and  V,  respectively. 
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The  extended  Kalman  filter  equations  are  similar  to  the  Kalman  equations,  as 
shown  in  the  following  equations. 

x(k+l|k)=f(x(k|k),k)  (19) 

z(k  +  l|k)  =  h(x(k  +  l|k),k)  (20) 

x(k+l|k+l)=x(k+l|k)+G(k+l)[z(k+l)-z(k+l|k] ,  (21) 

where:  x(k|k)  =  the  current  state  estimation, 

x(k  +  l|k)  =  the  predicted  state  estimation, 
z(k  +  l|k)  =  the  predicted  measurement 
z(k  + 1)  =  the  observed  measurement 
G(k  + 1)  =  the  Kalman  gain 

The  state  is  predicted  by  using  (19)  and  the  current  estimation  of  the  state.  The 
noise  is  in  z(k  +  1).  Equation  (20)  predicts  the  next  measurement  using  the  estimate  of 
the  predicted  state  and  the  measurement  equation  h( ).  The  final  Equation  (21)  corrects 
the  estimate  of  the  state  ba^  upon  the  innovation  and  the  Kalman  gains.  The  innovation 
is  the  difference  between  actual  meastirements  and  the  predicted  measurements. 

The  Kalman  gains  must  be  computed  on-line'  in  the  extended  Kalman  filter  rather 
than  off-line  as  in  the  ordinary  Kalman  filter.  Here  is  where  the  linear  approximation  of 
the  nonlinear  model  enters  into  play.  First,  the  Jacobians 
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P(k  +  l|k  + 1)  =  [l  -  G(k  +  l)H]P(k  +  ijk) ,  (26) 

where:  P(k|k)=  the  estimated  covariance  matrix  for  the  system  state  following  the  k* 
me^urement  update. 

P(k  +  l|k)  =  the  estimation  of  the  covariance  matrix  for  the  predicted  state  prior 
to  measurement  update. 

Note  that  Equations  (19)  through  (26)  specify  the  discrete  extended  Kalman  filter 
equations.  Initial  estimates  must  be  provided  for  the  state  and  the  covariance  matrix. 
These  estimates  are  usually  chosen  as 

x(0|0)  =  E[x(0)]-.  (27) 

P(OlO)  =  E  |{x(0)  -  x(0|0)  }^(0)  -  x(OlO)  y  ^  Cov[x(0)] .  (28) 
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It  must  be  noted  that  since  (22)  and  (23)  are  linearized  about  the  current  estimate 
of  the  state,  the  Kalman  gains  must  be  computed  at  the  same  time  observations  are  being 
made;  in  real  time. 

The  Kalman  filter  guarantees  optimal  performance,  stability,  and  convergence. 
The  extended  Kalman  filter  makes  no  such  promises.  It  is,  in  fact,  a  suboptimal  filter. 
The  linear  approximation  that  the  Jacobian  evaluations  represent  can  lead  to  a  divergent 
filter. 
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V.  SAT  SIMULATION  RESULTS  WITHOUT  DELAY 


A.  EXTENDED  KALMAN  FILTER  EQUATIONS 

The  extended  Kalman  filter  will  be  applied  to  the  SAT  problem.  The  algorithm 
will  be  tested  by  simulating  it  in  MATLAB.  The  following  assumptions  are  made: 

1 .  The  locations  of  the  UAV,  the  van,  and  the  cell  base  station  are  known 
precisely. 

2.  A  system  will  be  in  place  to  measure  the  phase  differences  in  the  SAT 
from  the  cell  and  the  phone  and  convert  these  into  TDOA 
measurements.  To  this  measurement  a  white  Gaussian  noise  will  be 
added. 

3.  The  phone  is  not  moving. 

4.  All  locations  are  referenced  to  a  geostationary  coordinate  system. 

5.  The  delay  introduced  by  the  phone  in  returning  the  SAT  will  be 
neglected. 

These  assumptions  will  help  to  simplify  the  problem.  The  delay  introduced  by  the  phone 
will  be  handled  later  as  it  is  nontrivial.  The  state  and  observation  equations  follow. 
Referring  to  Figure  3: 


x(k  +  l)  =  x(k)=  , 

yt  K 

(29) 

Rbp  =  -  xcy + (^  -  yc^ , 

(30) 

RpR  =-J(xt-xry  +(yt-yTy  , 

(31) 

d  =  ^(xc-xr)^+(yc-yr)^ , 

(32) 
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(33) 


TDOA  =  Zr  =-[Rbp  +  RpR  -d]+v, 

where:  xt,  yt  =  the  estimate  of  the  phone's  coordinates, 

xr,  yr  =  the  receiver  coordinates  (the  van  or  UAV), 
xc,  yc  =  the  cell  base  station  coordinates, 
c  =  the  speed  of  light, 

V  =  white  measurement  noise. 

Notice  that  since  the  phone  is  stationary,  the  state  transition  matrix  is  the  identity 
matrix.  The  nonlinearity  of  the  observation  equation  (30-33)  is  what  leads  to  the  use  of 
the  extended  Kalman  filter.  In  implementing  the  extended  Kalman  filter,  each  receiver 
will  be  treated  separately.  .The  TDOA  for  the  van  will  be  calculated  from  the  current  van 
position,  and  the  estimated  phone  position.  This  will  be  used  by  the  filter  to  update  the 
predicted  phone  position,  which  is  fed  into  another  Kalman  filter  handles  the  TDOA  for 
the  UAV.  Each  filter  works  off  the  results  of  the  prior  filter.  In  this  way,  the  two  stage 
filter  should  converge  estimates  upon  the  true  position  of  the  phone.  A  two  stage  filter 
algorithm  is  used,  rather  than  trying  to  combine  both  into  a  single  stage  filter. 

To  use  the  extended  Kalman  algorithm,  the  Jacobian  of  h(  )  must  be  foimd. 
Taking  the  derivatives  and  making  substitutions, 


hx=-!^ 

c 


xt-xc  xt-xr  xr-xc 
•+ - 


R 


BP 


"PR 


(34) 
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(35) 


1  yt-yc_^yt-yr  yr-yc 
c  _  R  BP  R  pp  d 


H  =  [hxhy].  (36) 

Equation  (36)  can  then  be  used  in  the  extended  Kalman  filter  equations 

P(k+llk)=P(k|k>Q,  (37) 

P(k+l|k+l)=[l-G(k+l)H(k)]p(k+l|k),  (38) 

G(k+l)=P(k+l|k)H(ky[H(k)P(k+l|k)H(ky+R]’'  .  (39) 

Note:  H(k)whouldbe  H(k+llk)ifthe  phone  was  not  stationary. 


For  stability,  the  Joseph  form  of  the  covariance  update  equation  will  be  used 
instead  of  Equation  (38). 

B  =  [l-G(k  +  l)H(k)],  (40) 

P(k+l|k+l)= BP(k+l|k)B^  +G(k+l)RG(k+l)  \  (41) 

The  use  of  the  Joseph  form  is  more  computationally  expensive,  but  is  less 
sensitive  to  round-off  errors  and  will  not  lead  to  negative  eigenvalues. 
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B.  TDOA  SIMULATIONS 


To  test  the  algorithm  the  following  scenario  was  used.  The  cell  base  station  is 
placed  at  the  coordinates  (10,10).  The  intelligence  van  is  placed  at  (10,0).  The  UAV  is 
free  to  fly  a  circular  pattern  that  will  be  centered  either  on  the  cell  center  or  the  van.  The 
phone  is  free  to  be  placed  inside  or  oUtside  the  pattern  of  the  UAV.  The  initial  error 
covariance  matrix  and  the  variance  of  plant  excitation  are  varied  until  the  filter  reaches  its 
best  convergence  characteristics.  Arbitrary  white  noise  has  been  added  to  the  TDOA 
measurements.  The  initial  estimation  of  the  phone  position  will  always  be  (8,8).  There  is 
no  particular  reason  for  this  and  the  initial  guess  is  not  plotted. 

For  each  scenario  three  plots  are  shown.  The  first  shows  a  time  vs.  state 
estimation  for  the  phone.  This  shows  the  convergence  characteristics  of  the  filter.  The 
second  plot  will  be  an  X-Y  plot  showing  the  trajectory  of  the  estimates  of  phone  position. 
This  plot  will  also  show  the  estimated  loci  of  position  for  every  fifth  observation.  This 
gives  an  idea  of  how  orthogonal  the  loci  are.  The  UAV  will  be  represented  by  an 
asterisk,  the  track  of  the  estimated  position  by  a  jagged  solid  line,  and  the  loci  of  position 
by  dashed  lines  or  ellipses.  The  ellipses  of  three  sigma  will  be  shown  as  a  solid  ellipse. 
The  third  plot  will  be  a  close  up  of  the  area  around  the  actual  phone  position  to 
investigate  how  accurate  the  filter  is. 

The  observations  were  first  taken  at  a  uniform  interval  on  the  UAV  trajectory. 
The  results  of  this  were  less  than  satisfying.  The  measurements  were  then  taken  at 
random  intervals  of  the  UAV  trajectory,  which  worked  much  better.  The  MATLAB 
programs  used  are  included  in  the  appendices. 
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1.  Scenario  One 


In  this  scenario,  the  UAV  is  allowed  to  fly  a  circle  around  the  cell  center  with  a 
radius  of  10  km.  The  phone  is  first  placed  at  (13,14).  The  best  convergence  occurred 
with 


Q= 


5  0 
0  5 


xl0^  P0= 


5  0 
0  5 


xl0‘®,  andR=l.694xlO 


.20 


(42) 


It  is  evident  in  Figure  7  that  the  filter  converges  quickly-within  8  observations. 
After  convergence,  the  filter  remains  stable. 

In  Figure  8  it  can  be  seen  that  the  trajectory  of  the  estimation  of  the  phone 
position  does  a  fairly  good  job  at  converging  on  the  phone  location.  It  does  show  some 
large  perturbances,  but  settles  down  on  the  phone  location  quickly.  Figure  9  shows  that 
the  filter  does  wander  around  the  location  of  the  phone  a  little.  The  final  two  three  sigma 
error  ellipses  are  shown  in  the  close-up.  These  show  that  there  is  a  probability  of  98% 
that  the  phone  is  in  this  ellipse.  The  final  ellipse  has  axis  of  dimension  180  m  by  40  m. 
This  is  very  acceptable. 
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Figure  7.  Observations  vs.  Phone  Location,  Scenario  One. 
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Estimate  of  phone 


Figure  8.  Trajectory  of  the  Phone  Estimation,  Scenario  One. 


(km) 


Figure  9.  Close-up  of  Trajectory,  Scenario  One. 
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2.  Scenario  Two 


The  next  test  keeps  everything  the  same  except  the  phone  position.  The  phone  is 
now  placed  at  (10,24),  outside  of  the  UAV  pattern  but  in-line  with  the  cell  center  and  the 
van.  The  results  are  shown  in  Figures  10,  1 1,  and  12. 

The  final  accuracy  of  the  filter  is  not  as  good  as  in  Scenario  One.  It  takes  the 
filter  just  over  20  observations  to  converge.  Figure  1 1  shows  how  the  ellipses  of  location 
intercept  each  other  around  the  phone  location.  The  intersection  is  not  as  orthogonal  as 
before,  which  would  cause  the  longer  convergence  time.  The  final  error  ellipse  is  slightly 
larger  than  before,  at  180  m  by  50  m.  The  filter  is  rather  well  behaved. 
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Phone  location  in  km 


Estimate  of  the  x  and  y  coordinates 


Figure  10.  Observations  vs.  Location,  Scenario  Two. 
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Figure  12.  Close-up  of  Trajectory,  Scenario  Two. 


3.  Scenario  Three 


As  a  final  test  the  phone  is  placed  at  (25,11).  The  filter  converges  within  5 
observations,  as  seen  in  Figure  13.  Figure  14  shows  that  the  intersection  of  the  ellipses  of 
location  are  fairly  orthogonal.  This  is  more  visible  in  the  close-up  of  Figure  15.  The 
major  axis  of  the  ellipse  is  aligned  parallel  to  the  location  ellipses,  which  is  expected. 
The  size  of  the  error  ellipse  is  30  m  by  350  m.  Again,  this  is  acceptable.  It  should  be 
noted  that  even  though  the  ellipse  has  such  a  large  major  axis,  the  actual  track  stays  well 
within  the  center  of  the  ellipse  and  over  the  phone's  position.  If  only  this  small  section  of 
the  track  is  inspected,  the  accuracy  is  quite  good.  Obviously  the  geometry  of  the  problem 
results  in  loci  of  position  that  are  not  as  orthogonal  as  the  two  other  situations. 
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Estimate  of  the  x  and  y  coordinates 


Figure  13.  Observations  vs.  Location,  Scenario  Three. 
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Figure  15.  Close-up  of  Trajectory,  Scenario  Three. 
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4.  Results 


The  algorithm  worked  fairly  well.  The  best  accuracy  occurred  when  the  phone 
was  placed  within  the  UAV  pattern.  The  worst  convergence  occurred  when  the  phone 
was  placed  outside  the  UAV  pattern  and  not  in-line  with  the  van  and  cell  base.  The 
ellipses  of  error  show  the  relative  nature  of  the  final  accuracy’s  for  the  filter. 

Observations  were  made  at  random  intervals  on  the  UAV  track,  which  helped  the 
filter  to  converge  both  faster  and  more  accurately.  Although  not  shown,  when 
observations  were  taken  at  fixed  intervals  of  its  flight  path,  the  filter  had  much  worse 
performance.  One  interesting  result  of  using  random  observations  is  that  each  simulation 
nm  on  the  same  initial  conditions  gave  different  results.  The  plots  shown  are 
representative  of  the  majority  of  these  runs. 
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VI.  SAT  SIMULATION  RESULTS  WITH  DELAY 
A.  EXTENDED  KALMAN  FILTER  EQUATIONS 

In  the  last  chapter  the  extended  Kalman  filter  was  applied  to  the  SAT  problem 
without  consideration  to  the  delay  encountered  at  the  phone.  Each  phone  has  some  delay 
associated  with  processing  signals.  This  is  usually  set  by  the  manufacturer.  If  the  delay 
was  known,  it  would  be  a  simple  matter  to  add  it  into  our  TDOA  equations.  It  must  be 
assumed  that  it  is  not  known,  and  as  an  unknown  variable  it  can  introduce  unacceptable 
error  into  our  filter.  Figure  16  shows  the  trajectory  that  results  when  using  the  filter 
developed  in  Chapter  V.  A  delay  of  4  microseconds  was  arbitrarily  chosen. 

The  position  of  the  phone  was  (13,14),  which  previously  yielded  good  results. 
The  filter  does  not  converge  upon  the  actual  phone  position,  but  is  offset  by  the  error 
introduced  by  the  delay.  Notice  that  the  y  coordinate  of  the  filter  converges  to  almost  15 
km.  The  close-up  view  in  Figure  17  further  shows  how  the  delay  affects  the  filter  result. 
It  does  not  converge  around  the  phone  as  before. 
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Estimate  of  the  x  and  y  coordinates 


Figure  16.  Trajectory  with  a  4  Microsecond  Delay. 
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To  compensate  for  the  unknown  delay,  the  filter  equations  must  be  modified.  The 
unknown  delay  is  now  considered  a  third  state  of  the  equation,  and  is  to  be  estimated. 
This  is  accomplished  rather  easily.  The  new  state  is 

xt(k) 

x(k+l)=x(k)=  yt(k)  .  (43) 

d(k) 

The  TDOA  equations  are  modified  by  adding  d  to  the  real  observations  and  d(k) 
to  the  estimated  TDOA.  The  Jacobian  is  also  modified  by 

H=[hxhyht],  (44) 

where  hx  and  hy  are  the  same  as  before. 

1.  Scenario  One 

The  initial  estimate  of  delay  will  be  5  microseconds.  The  initial  R  stays  the  same 
as  before.  The  challenge  with  this  algorithm  is  to  find  the  Q  and  PO  that  give  the  best 
convergence.  After  much  experimentation  the  following  was  used. 

'5x10*  0  0  1  [5x10*®  0  0 

Q=  0  5x10*  0  and  P0=  0  5x10'®  0  (45) 

0  0  2x10’“  0  0  2ixl0‘" 

The  results  of  this  filter  are  seen  in  Figures  18,  19,  and  20  for  a  phone  position  of 
(13,14)  and  an  actual  delay  of  4  microseconds.  Figure  18  shows  the  convergence  of  the 
states  within  10  observations.  In  Figure  19  it  can  be  seen  that  the  error  ellipsoids  start  out 
very  large  and  become  smaller  as  the  filter  converges.  The  close-up  in  Figure  20  shows 
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the  final  error  ellipsoid  has  the  dimensions  1.300  m  by  80  m.  This  is  larger  than  the 
results  shown  in  Chapter  V.  The  filter  converges  upon  the  phone's  position  rather  well. 
It  is  much  better  than  the  results  with  no  compensation  for  the  delay. 
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Estimate  of  the  x  and  y  coordinates 


Figure  18.  State  Plots  with  4  Microsecond  Delay,  Scenario  One. 
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2.  Scenario  Two 

Here  the  phone  is  placed  at  the  position  (10,24)  as  was  done  in  the  previous 
chapter.  The  results  are  seen  in  Figures  21  and  22. 

In  this  example,  the  delay  converges  to  25  microseconds  and  throws  off  the 
position  estimates.  The  y  coordinate  of  the  estimated  phone  position  converges  to  20  km, 
4  km  off  the  mark.  The  state  trajectories  in  Figure  22  show  Just  what  effect  the  error  in 
delay  estimation  has  on  the  position  estimation. 
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delay  in  microseconds  ,  Phone  location  in  km 


Estimate  of  the  x  and  y  coordinates 


Figure  21.  State  Plots  with  4  Mircosecond  Delay,  Scenario  Two. 
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3.  Scenario  Three 

The  phone  will  be  placed  at  (25,11)  which  gave  the  worst  results  for  the  filter 
discussed  in  Chapter  V.  The  results  are  seen  in  Figures  23,  24,  and  25.  The  filter 
converges  relatively  well,  coming  to  a  steady  state  within  10  observations.  The  delay 
state  is  still  not  converging  to  the  actual  state,  but  is  much  closer  than  before. 

The  major  axis  of  the  error  ellipsoid  is  large  at  over  2  km.  Note,  however,  that  the 
actual  track  converges  to  a  point  around  20  m  to  the  west  of  the  phone,  as  seen  in  Figure 
25.  If  one  was  looking  at  the  track  alone,  this  filter  would  give  acceptable  results. 
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Phone  location  in  km 


Estimate  of  the  x  and  y  coordinates 


Figure  23.  State  Plots  with  4  Microsecond  Delay,  Scenario  Three. 
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4.  Results 


The  above  plots  were  chosen  because  they  were  the  best  results  for  a  number  of 
runs  with  the  same  initial  conditions.  Each  run  gave  slightly  different  results.  This  filter 
gave  a  number  of  results  that  did  not  converge,  and  thus  the  filter  depends  upon  the 
random  measurements  for  its  stability.  While  the  number  of  results  that  did  not  converge 
was  not  great  the  filter  employed  in  Chapter  V  always  converged. 

By  treating  the  delay  as  a  state  to  be  estimated,  we  have  solved  the  problem  of 
most  of  the  offset  introduced  by  an  unknown  delay.  The  filter  did  not  always  converge 
on  the  actual  delay.  It  usually  got  close  enough  so  that  the  introduced  offset  error  was 
acceptable.  More  adjustment  of  the  Q  and  PO  matrices  may  make  the  delay  estimation 
always  converge  to  the  real  value. 
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VII.  FSK  BURST  TDOA  SIMULATION  RESULTS 


A.  EXTENDED  KALMAN  FILTER  EQUATIONS 

The  extended  Kalman  filter  is  applied  to  the  Burst  TDOA  problem.  The 
algorithm  will  be  tested  in  MATLAB  for  convergence  and  stability.  The  following 
assumptions  are  made: 

1 .  The  locations  of  the  UAV  and  the  van  are  known  precisely  (GPS  in 
each).  The  location  of  the  cell  antenna  is  not  needed. 

2.  A  signal  processing  system  will  be  in  place  to  filter  out  the  desired 
phone  registration  number  and  measure  the  arrival  time  to  compute  the 
TDOA  measurements.  To  this  measurement  a  white  Gaussian  noise 
measurement  will  be  added. 

3.  The  phone  is  not  moving. 

4.  All  locations  are  referenced  to  a  geostationary  coordinate  system. 

These  assumptions  help  to  simplify  the  problem.  The  estimated  state  equation  is 
the  same  as  before.  The  predicted  observation  equation  is  now 


z(k  +  l|k)  =  -(Rpv«-Rpu,v).  (46) 

c 

where:  Rp^3„  =  distance  between  estimated  phone  position  and  the  van,  and 
Rpu^=  distance  between  estimated  phone  position  and  the  UAV. 

Here  the  observed  TDOA  is  the  time  that  it  takes  the  coded  burst  from  the  phone 
to  reach  the  van  and  UAV.  The  Kalman  equations  are  the  same  as  those  employed 
before.  The  Jacobian  has  changed  slightly,  to 
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H(k)  =  -!- 

C 


f 


xt  -  xvan  xt  -  xuav 


R 


Pvan 


R 


Puav  ' 


yt  -  yvan  yt  -  yuav 


R 


Pvan 


R 


Puav 


(47) 


The  Joseph  form  of  the  covariance  update  is  again  used  to  minimize  round  off 
errors.  The  filter  algorithm  works  the  same  as  the  SAT  algorithm  did.  The  difference 
here  is  that  only  one  loci  of  estimated  position  can  be  calculated  for  each  TDOA 
measurement.  The  filter  only  has  one  stage,  which  processes  the  single  observation.  It 
then  loops  around  for  the  next  measurement.  It  is  noted  that  if  the  UAV  did  not  move, 
the  loci  of  estimated  position  would  not  move,  and  no  clear  interception  would  be  found. 
The  loci  would  also  not  be  very  orthogonal,  and  a  great  deal  of  error  would  be  present  if 
it  did  converge. 

Like  before  the  Q  and  PO  matrices  are  varied  until  the  filter  gives  the  best 
convergence.  The  van  and  UAV  locations  are  as  before.  The  phone  will  be  placed  in 
three  locations  as  before. 


B.  TDOA  SIMULATIONS 

In  this  case,  all  of  the  initial  conditions  will  be  the  same  as  before.  The  best 
convergence  occurred  with 


0 

lO’ 


and 


P0= 


'5x10^° 

0 


0 

5x10“ 


(48) 


The  simulation  is  allowed  to  run  for  40  observations.  The  observations  occur  at 
random  intervals  of  the  UAV  track.  This  not  only  gives  better  convergence,  but  simulates 
the  noisy  environment  and  difficult  detection  of  these  bursts.  The  three  scenarios  are  run 
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as  before  with  the  same  resulting  three  plots  shown.  The  lines  represent  the  loci  of 
estimated  position  in  the  vicinity  of  the  estimation.  These  will  show  how  orthogonal  the 
measurements  are. 

1.  Scenario  One 

The  phone  is  placed  inside  the  UAV  pattern  at  (13,14).  The  results  are  seen  in 
Figures  26,  27  and  28. 

This  filter  acted  very  stable.  It  is  seen  firom  Figure  26  that  the  filter  converges 
within  30  observations.  The  trajectory  shown  in  Figure  27  shows  a  no  nonsense  track 
and  heads  directly  for  the  phone's  position.  This  is  seen  clearer  in  Figure  28  as  the  track 
converges  directly  on  the  phone's  location.  Although  this  algorithm  is  slower  to 
converge,  it  appears  to  be  more  accurate  than  the  SAT  algorithm. 
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Number  of  obsaivations 

Figure  26.  Burst  TDOA  Problem,  Scenario  One. 
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Figure  28.  Close-up  of  Trajectory,  Scenario  One. 
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2.  Scenario  Two 

In  this  scenario  the  phone  is  placed  at  (10,24),  outside  of  the  UAV  pattern,  but 
opposite  the  van.  The  results  are  as  shown  in  Figures  29,  30,  and  31.  Figure  29  shows 
that  the  convergence  of  this  situation  is  much  slower  than  before.  It  can  be  seen  that  the 
filter  in  fact  has  not  fully  converged  by  the  60^**  observation.  Figure  30  shows  the  track 
heading  in  the  general  direction  of  the  phone,  getting  ever  closer.  The  close-up  in  Figure 
31  shows  that  the  track  approaches  to  within  10  m  of  the  phone's  position.  This  is  very 
accurate. 

The  filter  seems  to  try  to  converge  but  never  makes  it.  This  phone  position  also 
gave  a  great  number  of  results  which  did  not  converge  at  all. 
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Estimate  of  the  x  and  y  coordinates 


Figure  29.  Burst  TDOA  Problem,  Scenario  Two. 
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Estimate  of  phone 


Figure  30.  Trajectory  of  Burst  TDOA  Problem,  Scenario  Two. 
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Figure  31.  Close-up  of  Trajectory,  Scenario  Two. 
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3.  Scenario  Three 


The  final  test  places  the  phone  outside  of  the  UAV  pattern  at  (25,1 1).  The  results 
are  shown  in  Figures  32,  33,  and  34.  In  general,  the  new  position  gives  the  same  results 
as  found  in  Scenario  Two.  The  change  does  not  have  much  of  an  effect  upon  the  filter. 

This  filter  never  converges  around  the  exact  phone  location.  It  displays  the  same 
single  mindedness  by  constantly  getting  closer  to  the  phone. 
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Estimate  of  the  x  and  y  coordinates 


Figure  32.  Burst  TDOA  Problem,  Scenario  Three. 
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Figure  34.  Close-up  of  Trajectory,  Scenario  Three. 
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4.  Results 


The  algorithm  worked  quite  well  when  the  phone  was  placed  inside  the  UAV 
track  (Scenario  One).  The  random  nature  of  the  observations  again  caused  differing 
results  with  the  same  initial  conditions.  This  algorithm  has  the  advantage  that  the  cell  site 
does  not  have  to  be  known,  and  only  one  frequency  need  be  monitored.  The  results  are 
close  to  the  actual  phone  position.  In  fact,  when  the  phone  was  placed  inside  the  UAV's 
flight  circle,  the  results  were  very  accurate.  This  method  favors  having  the  UAV  fly 
around  the  suspected  position. 
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VIII.  COMBINED  SAT  AND  BURST  FILTER 


A.  ALGORITHM  CHANGES 

In  this  chapter  the  two  methods  of  finding  the  position  of  the  phone  will  be 
combined  into  one  filter.  This  assumes  that  both  the  SAT  information  and  the  Burst 
information  are  available,  and  that  the  location  of  the  cell  station  is  known.  This  method 
is  easily  implemented  by  taking  the  last  two  filters  and  combining  them  in  series. 
Basically,  in  one  loop  of  the  program,  the  filter  will  first  modify  its  position  estimate 
based  upon  the  SAT  information.  It  will  then  feed  this  estimate  to  the  filter  produced  for 
the  Burst  filter.  The  algorithm  then  loops  back  to  the  beginning  using  the  Burst  updated 
estimation.  This  is  certainly  the  easiest  implementation  of  this  idea.  The  program  could 
also  have  been  modified  to  use  both  measurements  of  TO  A  at  the  same  time.  To  keep 
things  simple,  it  is  assumed  that  the  SAT  delay  is  zero  for  the  SAT  calculations.  This  has 
been  done  to  test  if  these  two  algorithms  can  work  together  on  a  basic  level. 

B.  SIMULATIONS 

The  simulations  will  follow  the  same  pattern  used  in  the  previous  chapters.  The 
same  locations  are  used  so  that  the  results  may  be  compared  to  prior  results.  Random 
measurements  will  again  be  assumed  on  the  track  of  the  UAV.  This  is  not  so  unrealistic, 
since  in  a  real  scenario  the  receivers  may  be  able  to  discern  only  occasional  pieces  of 
information  firom  the  cellular  phone.  As  the  phone  travels  through  changing  landscape, 
each  receiver  may  temporarily  find  itself  in  a  blind  spot  of  the  signal.  This  algorithm  will 
test  how  well  one  can  predict  the  phone  position  when  the  SAT  and  Burst  information  are 
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spotty.  Hopefully,  the  combined  information  will  be  enough  to  let  the  combined  filter 
converge. 

1.  Scenario  One 

The  phone  is  placed  at  (13,14)  to  test  how  the  algorithm  reacts  to  having  the 
phone  inside  of  the  UAV  pattern.  The  results  can  be  seen  in  Figures  35,  36,  and  37.  In 
Figure  35  it  is  seen  that  the  filter  converges  within  4  observations.  This  is  very  quick  in 
relation  to  the  other  filters.  Figure  36  shows  how  the  filter  heads  straight  for  the  phone's 
location.  The  close-up  in  Figure  37  gives  a  better  idea  of  what  is  happening  around  the 
point  of  convergence.  The  three  sigma  ellipses  are  160  m  by  100  m.  This  is  not 
significantly  smaller  than  the  other  results,  but  when  the  plot  is  examined  it  is  evident  that 
the  filter  stays  within  20  m  of  the  phone.  This  sort  of  accuracy  is  what  is  needed  for 
practical  application. 
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Estimate  of  the  x  and  y  coordinates 


Figure  35.  Burst  TDOA  Problem  with  Combined  Filter,  Scenario  One. 
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2.  Scenario  Two 


As  in  the  previous  chapters,  the  phone  is  placed  at  (10,24)  to  test  how  the 
algorithm  handles  the  situation  with  the  phone  outside  of  the  UAV  track  and  in-line  with 
the  cell  center  and  the  van. 

Figure  38  shows  the  filter  converging  within  15  observations.  It  is  seen  in  Figure 
39  that  the  filter  makes  the  largest  jumps  within  the  first  4  observations.  Figure  40,  the 
close-up  of  the  phone,  shows  the  three  sigma  ellipse  to  be  200  m  by  70  m.  The  final 
convergence  does  not  leave  the  phone's  position  by  more  than  40  m  once  it  has  settled 
down. 
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Estimate  of  the  x  and  y  coordinates 


Figure  38.  Burst  TDOA  Problem  with  Combined  Filter,  Scenario  Two. 
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3.  Scenario  Three 


In  this  scenario  the  phone  is  placed  at  (25,11)  to  test  the  response  of  having  the 
phone  outside  the  UAV's  track  and  not  in  line  with  the  van  and  cell  station. 

Figure  41  shows  that  this  filter  converges  quickly  in  this  scenario,  at  less  than  10 
observations.  The  track  of  the  filter  shown  in  Figure  42  shows  that  the  filter  has  some 
trouble  in  the  first  couple  of  observations,  estimating  the  position  to  be  as  far  as  65  km 
away  from  the  actual  position.  What  is  also  interesting  is  that  the  estimates  stay  outside 
of  the  UAV  track.  The  close-up  shown  in  Figure  43  shows  the  three  sigma  ellipse  having 
the  dimensions  of  160  m  by  120  m.  Once  converged  the  algorithm  keeps  its  estimates 
within  40  m  of  the  phone's  actual  position. 
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Estimate  of  the  x  and  y  coordinates 


Figure  41 .  Burst  TDOA  Problem  with  Combined  Filter,  Scenario  Three. 
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Close  up  of  position  of  the  Phone 
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C.  RESULTS 


The  biggest  benefit  of  combining  the  SAT  and  Burst  filters  is  not  shown  in  the 
plotted  figures,  as  the  combined  filter  converged  every  time  it  was  run.  The  stability  of 
this  filter  far  surpasses  the  stability  of  its  two  components  working  separately.  Out  of  the 
50  runs  made  with  this  algorithm,  only  one  failed  to  converge.  Remember  that  the  point 
of  the  reading  is  taken  at  a  random  spot  on  the  UAV  track  so  that  each  run  is  running  off 
of  essentially  unique  data  each  time.  The  robustness  of  this  algorithm  is  of  great  value. 
The  accuracy  and  time  to  converge  is  also  as  good  if  not  better  than  the  other  two  filters. 
It  converges  as  fast  as  the  SAT  filter  and  has  as  good  accuracy  as  either  filter. 

It  should  be  noted  that  both  an  SAT  and  Burst  data  reading  was  taken  at  each 
random  measurement  point.  In  the  future,  each  point  may  be  limited  as  to  which  signal 
can  be  received.  The  filter  would  have  to  be  modified  to  operate  only  when  data  was 
received,  ignoring  the  missing  data  in  the  loop.  The  delay  estimation  would  also  have  to 
be  included  in  the  next  iteration  of  this  program,  as  the  nonstandard  delay  of  the  SAT 
signal  can  affect  accuracy. 
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IX.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

All  of  the  extended  Kalman  filter  algorithms  developed  here  perform  well.  The 
position  of  a  cellular  phone  can  be  estimated  to  a  great  degree  of  accuracy  by  measuring 
the  TDOA  for  certain  cellular  emissions.  The  use  of  a  base  station  and  a  UAV  has  shown 
to  be  highly  beneficial  to  these  algorithms.  The  movement  of  the  UAV  causes  the 
shifting  of  the  loci  of  positions  for  the  filters.  This  shifting  of  the  loci  allows  the  Kalman 
filter  to  estimate  the  intersection  and  thus  the  phone  position. 

The  filters  developed  are  highly  dependent  upon  the  type  of  emission  being 
exploited.  Each  type  of  emission  has  strengths  and  weaknesses  with  regard  to  locating 
phones.  The  SAT  filter  is  very  accurate,  and  has  a  good  convergence  time.  The  SAT 
signal  is  continuous  and  easy  to  intercept  once  the  channels  being  used  are  known.  The 
major  detraction  is  that  the  phone  has  an  inherent  delay  which  adds  an  error  constant. 
This  thesis  has  shown  that  it  is  possible  to  predict  this  delay  by  using  a  three  state 
Kalman  filter.  This  filter  worked  fairly  well,  but  was  the  most  unstable  filter  of  the 
group.  Many  times  it  would  diverge  on  its  estimates  or  converge  to  the  wrong  delay 
value,  which  does  not  solve  the  delay  problem.  The  filter  does  show  that  it  can  work 
well,  though,  and  has  the  potential  to  become  a  very  stable  and  accurate  filter,  given 
further  development. 

The  Burst  Kalman  filter  relies  on  the  reception  of  control  codes  from  the  phone. 
This  reception  and  isolation  of  the  signal  is  a  thesis  problem  in  itself.  Assuming  that  this 
information  can  be  obtained,  the  filter  showed  that  the  concept  of  measuring  only  the 
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TDOA  does  work  well.  This  filter  showed  good  stability,  but  was  slow  to  converge  and 
often  never  actually  reached  the  phone  position. 

When  the  two  filters  were  combined,  ignoring  the  delay  of  the  SAT  problem,  the 
resulting  filter  had  the  best  features  of  both  of  its  components.  The  SAT  Burst  filter  was 
accurate,  quick  to  converge,  and  very  stable.  The  results  of  each  filter  enhances  the 
performance  of  the  next.  This  combination  shows  the  most  promise  since  one  or  the 
other  signal  or  both  may  be  available  at  any  one  point  in  time. 

These  algorithms  have  the  advantage  that  the  receivers  are  very  simple  and 
omnidirectional.  The  UAV  is  perfect  for  such  use,  since  the  hardware  required  by  such  a 
system  would  easily  fit.  (Though  a  small  plane  or  boat  could  also  be  used).  The 
algorithms  are  not  too  complex  and  would  not  overpower  a  small  mobile  computer. 

B.  RECOMMENDATIONS 

This  thesis  is  just  the  beginning  of  an  application  which  has  great  potential. 
Although  only  one  UAV  pattern  was  investigated  here,  different  patterns  may  have 
differing  efficiencies.  A  development  of  these  patterns  and  their  uses  would  have  merit. 
The  study  of  how  well  these  filters  track  moving  targets  would  be  a  natural  progression. 
The  a  priori  estimates  have  a  great  effect  upon  the  filter's  behavior.  A  study  of  how  these 
estimates  affect  stability,  convergence,  and  accuracy  would  be  worthwhile. 

It  would  be  beneficial  to  develop  an  algorithm  that  would  be  able  to  set  the  a 
priori  estimates  based  upon  expected  conditions  and  the  first  observation.  A  neural  net 
might  be  trained  to  do  this  so  that  the  filter  gives  optimum  results. 
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Two  receivers  could  have  been  used  on  each  vehicle.  This  might  increase 
accuracy.  Two  UAV's  might  also  be  used,  and  their  interaction  would  be  an  interesting 
area  of  investigation. 

The  combined  SAT  Burst  filter  should  be  modified  to  account  for  SAT  delay  at 
the  cellular  phone  and  to  operate  on  available  data.  Some  data  may  be  suppressed  or 
erroneous.  The  algorithm  should  be  able  to  distinguish  which  is  useable  and  apply  it  to 
the  Kalman  filter. 

When  all  is  said  and  done  this  idea  and  technique  show  merits  for  the  tracking  of 
cellular  phones.  The  widespread  use  of  these  phones  makes  this  tracking  even  more 
attractive. 
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APPENDIX  A.  ERROR  ELLIPSES 

In  Ref.  3,  a  theory  is  developed  that  derives  the  error  ellipsoids  from  a  given 
covariance  matrix  of  estimation  error.  This  is  P  in  this  work.  The  source  states  that  the 
eigenvalues  of  P  are  a,,a2,...,a„  and  the  corresponding  eigenvectors  are 

.  The  quadratic  form  of  the  hyperellipsoid  is  developed  as 

a,  ttj  a„ 

This  forms  an  ellipsoid  in  n  dimensional  space,  where  n  corresponds  to  the 
number  of  states  our  filter  has  and  y*  is  the  linear  coefficient  of  y^*^  and  so  on.  If  n  is 
greater  than  two  hyperellipsoids  are  formed.  These  hyperellipsoids  have  surfaces  which 
have  equal  probability  densities.  Thus,  for  a  certain  f  the  probability  that  a  point  will  lie 
inside  the  ellipsoid  can  be  computed. 

The  probabilities  have  been  calculated  for  a  few  values  of  dimension,  n  and  ^ ; 


1 

n 

1 

2 

3 

1 

0.683 

0.955 

0.997 

2 

0.394 

0.862 

0.989 

3 

0.2 

0.739 

0.971 

Table  Al.  Probability  ol 

'  n  Dimension  1 

Ellipsoids. 
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Since  it  is  desired  to  show  the  accuracy  potential  of  the  filters  and  the  dimension 
is  n=2  for  plotting  purposes,  1  =  3  or  a  three  sigma  error  ellipse  is  used.  This  final 
ellipse  yields  an  ellipse  for  which  there  is  a  99.7%  probability  that  the  desired  value  is 
inside. 

The  MATLAB  program  which  does  this  uses  the  MATLAB  "eig"  function  to 
calculate  the  eigenvectors  and  eigenvalues  of  the  covariance  matrix.  It  then  forms  the 
ellipse  from  Equation  lA  and  converts  this  into  a  polar,  which  is  easier  to  get  locations 
for  plotting. 
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APPENDIX  B.  MATLAB  PROGRAMS 


TDOA  -  THE  BURST  FILTER 
%  Michael  P.  Fallon 

%  TDOA;  This  program  simulates  the  TDOA  burst  estimation  problem. 
%  The  van  and  UAV  are  assumed  to  have  single  antennas  in 
%  this  scenario.  The  UAV  is  allowed  to  fly  a  racetrack  pattern 
%  The  reception  of  the  identification  signals  from  the  phone  are 
%  assumed  to  occur  at  random  intervals  on  the  UAV  track. 

% 

% 


clear 

c=3.00e+08;  %  speed  of  light 

R1  =1 .694e-20;  %  The  covariance  of  measurement  error 

Q=1 0e+04*eye(2):  %  Variance  of  Plant  excitation  noise 
P0=5e+20*eye(2);  %  Initial  error  covariance 

%  Choose  to  calculate  the  error  covariance  of  3  sigma 
C=3.0; 


tf=60: 

t=(0:tf); 


%  the  time  vector 
len=length(t); 
l=ones(1  ,len); 

%  Initialize  storage  for  Kalman  gains  and  the  states 
g=zeros(2,len); 

X=zeros(2,len); 

%  Emitter  location  input  by  the  user. 

xt=input(’  Enter  the  x  coordinate  of  the  phone  in  km..’); 

yt=input(’  Enter  the  y  coordinate  of  the  phone  in  km..’); 

xt=xt*1000; 

yt=yt*1000; 

%  A  priori  estimate  of  the  location  of  the  target 
X0=[1 0000; 10000]; 

%  Location  of  the  van  and  UAV 
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xvan=10000*l+0*t; 

yvan=0*l+0*t; 

%  The  UAV  will  move  in  a  circular  pattern  centered  on  the  cell. 

dtheta=pi*rand(1  ,len)/5: 

theta(1)=-pi/2; 

for  i=2:len 

theta(i)=theta(i-1  )+dtheta(i); 
end 

%  Set  the  center  of  rotation  and  the  radius  for  the  UAV. 
Cx=10000; 

Cy=10000: 

r=8000; 

xuav=Cx*l+(r*l).*cos(theta); 

yuav=Cy*l+(r*l).*sin(theta); 

%  Calculate  the  TDOA  observations  for  the  two  receivers 
Z=(1/c)*(sqrt((xt-xvan).^2+(yt-yvan).'^2)-... 
sqrt((xt-xuav).'^2+(yt-yuav).'^2)): 

%  Inject  zero  mean  noise  into  the  observations 
n=sqrt(R1  )*randn(1  ,len); 

Zn=Z+n; 

%  Initialize  the  Kalman  Filter 

X(: ,  1  )=X0;  %  A  priori  estimate  of  the  emitter  location 

PK=PO:  %  Initial  error  covariance  matrix 

%  calculate  the  filtered  estimates  of  the  emitter  location, 
for  k=1:len-1, 

Rvan=sqrt((X(1  ,k)-xvan(k)).A2+(X(2,k)-yvan(k)).A2); 
Ruav=sqrt((X(1,k)-xuav(k)).^2+(X(2,k)-yuav(k)).'^2); 

hx=(1/c)*((X(1  ,k)-xvan(k))./Rvan-(X(1  ,k)-xuav(k))./Ruav); 
hy=(1/c)*((X(2,k)-yvan(k))./Rvan-(X(2,k)-Vuav(k))./Ruav); 

HKi=[hx  hy]; 

%  calculate  the  estimate  of  TDOA  based  upon  the  estimate  of 
%  phone  location  and  location  of  van  and  UAV 
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Zhat=(1/c)*(sqrt((X(1,k)-xvan(k))^2+(X(2,k)-yvan(k))A2)-... 

sqrt((X(1,k)-xuav(k))^2+(X(2,k)-yuav(k))^2)); 

Zh(k)=Zhat; 

%  The  new  error  covariance  is  the  same  as  the  old  plus  Q 
PK=PK+Q: 

%  Calculate  the  Kalman  gains 
GK=PK*HK’*inv(HK*PK*HK’+R1 ); 

%  Calculate  the  updated  error  covariance  matrix 
%PK=(eye(2)-GK*HK)*PK: 

%  Joseph  form  of  covariance  update. 
PK=(eye(2)-GK*HK)*PK*(eye(2)-GK*HK)’+GK*RrGK’; 

%  Calculate  the  smoothed  estimate  of  the  bearing  to  the  phone 
X(:,k)=X(:,k)+GK*(Zn(k)-Zhat); 

X(:,k+1)=X(:,k); 


%  save  the  error  covariance  matrices 

if  k==1 

P=PK; 

else 

P=[P',PKl 

end; 

%  save  the  Kalman  gains 
g(:.k)=GK; 

end; 


%  Plot  the  output 

%  Plot  the  estimates  of  the  phone  location. 

Figure(1) 

cig; 

plot(t,X(1  ,:)/1 000,t,X(2,:)/1 000) 

titleC Estimate  of  the  x  and  y  coordinates’); 

xlabel(’Number  of  observations’) 

ylabeli’ Phone  location  in  km’) 

grid 

pause 
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%  Steady  state  estimation  of  emitter  location 
xtavg=mean(X(1  ,k-10:k)); 
ytavg=mean(X(2,k-1 0:k)); 

%  Plot  the  grid 
Figure(2) 

plot(xt/1000,yt/1000,’*’) 
hold  on 

plot(xvan(k)/1 000,yvan(k)/1  OOO/r*’) 
plot(X0(l  )/1 000,X0(2)/1 000,’+’) 

%  Plot  the  error  ellipses, 
for  k=1:len-1; 

piot(xuav(k)/1 000,yuav(k)/1 000, ’g*’); 

P1=[P(2*k-1,1)  P(2*k-1,2);P(2*k,1)  P(2*k,2)] 
[xg,yg,Emax]=ellipa(P1  ,C,X(1  ,k),X(2,k)); 

%plot(xg/1 000,yg/1 000,’b-’); 
end 

plot(X(1  ,:)./1 000,X(2,:)./1 000, ’r’); 

axis(’equar) 

hold  off 

title(’Estimate  of  phone’) 
xlabel(’(km)’);ylabel(’(km)’); 

pause 

%  Plot  close  up  of  the  estimate  of  the  emitter  location 

Figure(3) 

cig 

axis([(xt-2000)/1000  (xt+2000)/1000  (yt-2000)/1000 ... 

(yt+2000)/1000]); 

hold  on 

plot(xt/1000,yt/1000,’*’) 
plot(X(1  ,:)./1 000,X(2,:)./1 000); 


k=len-2; 

p=P(2*k-1:2*k,:); 

[xg,yg,Emax]=ellipa(p,C,xtavg,ytavg); 
%plot(xg/1 000,yg/1 000,’b-’); 

[xg,yg,Emax]=ellipa(P1,C,xtavg,ytavg); 
%plot(xg/1 000,yg/1 000,’b-’); 
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hold  off 

title(’Close  up  of  position  of  the  Phone’) 
xlabel(’(km)’):ylabel(’(km)’); 

TDOA2-SAT  FILTER  WITHOUT  DELAY 


%  Michael  P.  Fallon 

%  TDOA2  This  program  simulates  the  TDOA  estimation  problem. 
%  In  this  scenario  we  assume  that  the  van  and  UAV  are  using 
%  the  SAT  to  calculate  the  ellipses  of  location. 

clear 

c=3.00e+08: 

%  speed  of  light 

R1  =1 .694e-1 5;  %  The  covariance  of  measurement  error 

Q=5e+02*eye(2);  %  Variance  of  Plant  excitation  noise 

P0=5e+1 0*eye(2):  %  Initial  error  covariance 

%  Choose  to  calculate  the  error  covariance  of  3  sigma 
C=3; 


tf=50; 

t=(0:tf); 


%  the  time  vector 
len=length(t); 
l=ones(1  ,len); 

%  Initialize  storage  for  Kalman  gains  and  the  states 
g=zeros(2,len); 

X=zeros(2,len); 

%  Emitter  location 

xt=input(’  Enter  the  x  coordinate  of  the  phone  in  km..’); 
yt=input(’  Enter  the  y  coordinate  of  the  phone  in  km..’); 
d=input(’  Enter  the  phone  delay  in  microseconds..’); 
d=d*1e-6; 

xt=xt*1000; 

yt=yt*1000; 
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%  A  priori  estimate  of  the  location  of  the  target 
X0=[8000:8000]; 

%  Location  of  the  cell  site  antennae 

xc=1 0000*1; 
yc=1 0000*1; 

%  Location  of  the  van  and  UAV 
xvan=1 0000*l+0*t; 
yvan=0*l+0*t; 

dtheta=pi*rand(1  ,len)/2; 

%dtheta=pi*.1; 

theta(1)=-pi/2; 

for  i=2:len 

theta(i)=theta{i-1  )+dtheta(i); 
end 

Cx=10000; 

Cy=10000; 

r=10000; 

xuav=Cx*l+(r*l).*cos(theta); 

yuav=Cy*l+(r*l).*sin(theta); 

%xuav=1 5*1+1  *t; 

%yuav=0*l; 

%  Calculate  the  TDOA  observations  for  the  two  receivers 

Z1=(1/c)*(sqrt((xt-xc).^2+(yt-yc).'^2)... 

+sqrt((xt-xvan).^2+(yt-yvan).^2)... 

-sqrt((xvan-xc).''2+(yvan-yc).'^2))+d; 

Z2=(1  /c)*(sqrt((xt-xc).^2+(yt-yc).^2) ... 

+sqrt((xt-xuav)  .^2+(yt-yuav).^2). . . 

-sqrt((xuav-xc)  .'^2+(yuav-yc)  .^2))+d; 

%  Inject  zero  mean  noise  into  the  observations 

n  1  =sqrt(R1  )*randn(1  ,len); 
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n2=sqrt(R1  )*randn(1  Jen); 


Z1n=Z1+n1; 

Z2n=Z2+n2: 

%  Initialize  the  Kalman  Filter 

X(:,1  )=X0;  %  A  priori  estimate  of  the  emitter  location 

PK=PO;  %  Initial  error  covariance  matrix 

%  calculate  the  filtered  estimates  of  the  emitter  location 
for  k=1:len-1, 

%  perform  the  first  calculation  for  the  measurement  at  the  van. 

Rtc=sqrt((X(1,k)-xc(k)).A2+(X(2,k)-yc(k)).A2): 

Rtr=sqrt((X(1  ,k)-xvan(k)).'^2+(X(2,k)-yvan(k)).^2); 
Rrc=sqrt((xvan(k)-xc(k)).^2+(yvan(k)-yc(k)).''2); 

h1  x=(1/c)*((X(1  ,k)-xc(k))./Rtc+(X(1  ,k)-xvan(k))./Rtr-... 
(xvan(k)-xc(k))./Rrc); 

h1y=(1/c)*((X(2,k)-yc(k))./Rtc-(X(2,k)-yvan(k))./Rtr-... 

(yvan(k)-yc(k))./Rrc); 

HK=[h1xh1y]; 

%  calculate  the  estimate  of  TDOA  based  upon  the  estimate  of 
%  phone  location  and  location  of  van  and  UAV 

Zhatl  =(1/c)*(sqrt((X(1  ,k)-xc(k)).'^2+(X(2,k)-yc(k)).^2)... 
+sqrt((X(1,k)-xvan(k)).''2+(X(2,k)-yvan(k)).^2)... 
-sqrt((xvan(k)-xc(k)).^2+(yvan(k)-yc(k)).'^2)); 

%  The  new  error  covariance  is  the  same  as  the  old  plus  Q 
PK=PK+Q; 

%  Calculate  the  Kalman  gains 
GK=PK*HK’*inv(HK*PK*HK’+R1 ); 

%  Calculate  the  updated  error  covariance  matrix 
%PK=(eye(2)-GK*HK)*PK; 

%  Joseph  form  of  covariance  update. 
PK=(eye(2)-GK*HK)*PK*(eye(2)-GK*HK)’+GK*RrGK’; 
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%  Calculate  the  smoothed  estimate  of  the  bearing  to  the  phone 
X(:,k)=X(:,k)+GK*(Z1  n(k)-Zhat1 ); 

X(;,k)=X(:,k): 

%  perform  the  calculations  for  the  measurement  at  the  UAV 

Rtc=sqrt((X(1,k)-xc(k))A2+(X(2.k)-yc(k))A2); 

Rtr=sqrt((X(1  ,k)-xuav(k)).'^2+(X(2,k)-yuav(k)).''2): 
Rrc=sqrt((xuav(k)-xc(k)).^2+(yuav(k)-yc(k)).^2); 

h2x=(1/c)*((X(1  ,k)-xc(k))./Rtc+(X(1  ,k)-xuav(k))./Rtr-... 
(xuav(k)-xc(k))./Rrc); 

h2y=(1/c)*((X(2,k)-yc(k))./Rtc-(X(2,k)-yuav(k))./Rtr-... 

(yuav(k)-yc(k))./Rrc); 

HK=[h2x  h2y]; 

%  calculate  the  estimate  of  TDOA  based  upon  the  estimate  of 
%  phone  location  and  location  of  van  and  UAV 

Zhat2=(1  /c)*(sqrt((X(1  ,k)-xc(k))  A2+(X(2,k)-yc(k))  .^2)... 

+sqrt((X(  1  ,k)-xuav(k))  .^2+(X(2,  k)-yuav(k)).^2). . . 
-sqrt((xuav(k)-xc(k))A2+(yuav(k)-yc(k)).''2)); 

%  The  new  error  covariance  is  the  same  as  the  old  plus  Q 
PK=PK+Q; 

%  Calculate  the  Kalman  gains 
GK=PK*HK’*inv(HK*PK*HK’+R1 ); 

%  Calculate  the  updated  error  covariance  matrix 
%PK=(eye(2)-GK*HK)*PK; 

%  Joseph  form  of  covariance  update. 
PK=(eye(2)-GK*HK)*PK*(eye(2)-GK*HK)’+GK*R1*GK’; 

%  Calculate  the  smoothed  estimate  of  the  bearing  to  the  phone 
X(:,k)=X(:,k)+GK*(Z2n(k)-Zhat2); 

X(:,k+1)=X(:,k); 


if  k==1 

Zhat=[Zhat1;Zhat2]; 

else 

Zhat=[Zhat;Zhat1  ;Zhat2]; 
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end 


%  save  the  error  covariance  matrices 

if  k==1 

P=PK: 

else 

P=[P:PK]: 

end; 

%  save  the  Kalman  gains 
g(:.k)=GK; 

end; 


%  Plot  the  output 

Figure{1) 

cig; 

plot(t,X(1  ,:)/1 000,t,X(2,:)/1 000) 

title(’ Estimate  of  the  x  and  y  coordinates’); 

xlabel(’Number  of  observations’) 

ylabel(’Phone  location  in  km’) 

grid 

pause 

%  Steady  state  estimation  of  emitter  location 
xtavg=mean(X(1  ,len-1 0:len)); 
ytavg=mean(X(2,len-1 0:len)); 


%  Plot  the  grid 
Figure(2) 

plot(xt/1000,yt/1000,’*’) 
hold  on 

plot(xvan(k)/1 000,yvan(k)/1 000, ’r*’) 
plot(X0(1  )/1 000,X0(2)/1 000,’+’) 


for  k=1  :len-1 ; 

plot(xuav(k)/1 000,yuav(k)/1 000,’g*’); 
P1=:[P(2*k-1,1)  P(2*k-1,2);P(2*k,1)  P(2*k,2)]; 
[xg,yg,Emax]=ellipa(P1  ,C,X(1  ,k),X(2,k)); 
plot(xg/1000,yg/1000,’b-’); 
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if  rem(k,10)==0 

[evan]=eliip1  (xc(k),yc(k),xvan(k),yvan(k),Zhat(2*k-1  ),c); 

[euav]=ellip1(xc(k),yc(k),xuav(k),yuav(k),Zhat(2*k),c); 

plot(evan(1,:)/1000,evan(2,:)/1000,’g:’): 

plot(euav(1  ,:)/1 000,euav(2,:)/1 000, ’y:’); 

end 

end 

plot(X(1  000,X(2,:)./1 000, ’r’); 

axis(’equar) 

hold  off 

title(’Estimate  of  phone’) 

xlabel(’(km)’):ylabel(’(km)’); 

pause 

%  Plot  close  up  of  the  estimate  of  the  emitter  location 

Figure(3) 

cig 

axis([xt/1 000-1  xt/1 000+2  yt/1 000-1  yt/1 000+2]) 
hold  on 

plot(xt/1000,yt/1000,’*’) 

plot(X(1,:)./1000,X(2,:)./1000); 

k=len-2; 

p=P(2*k-1:2*k,:): 

[xg,yg,Emax]=elllpa(p,C,X(1,len-2),X(2,len-2)); 

plot(xg/1000,yg/1000,’b-’); 

[xg,yg,Emax]=ellipa(P1  ,C,X(1  ,len-1  ),X(2,len-1 )); 
plot(xg/1 000,yg/1 000,’b-’); 
plot(xtavg/1 000,ytavg/1 000,’r+’): 

title(’Close  up  of  position  of  the  Phone’) 
xlabel(’(km)’);ylabel(’(km)’); 

hold  off 


TDOA2A  -SAT  FILTER  WITH  DELAY 
%  Michael  P.  Fallon 

%  TDOA2A  This  program  simulates  the  TDOA  estimation  problem. 
%  In  this  scenario  we  assume  that  the  van  and  UAV  are  using 
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%  the  SAT  frequency  to  calculate  the  ellipses  of  location 
clear 

c=3.00e+08:  %  speed  of  light 

R1  =1 .694e-1 5;  %  The  covariance  of  measurement  error 

Q=5e+04*eye(3):  %  Variance  of  Plant  excitation  noise 

Q(3,3)=2e-25: 

P0=5e+8*eye(3);  %  Initial  error  covariance 

P0(3,3)=2.5e-11: 

%  Choose  to  calculate  the  error  covariance  of  3  sigma 
C=3..0: 


tf=30: 

t=(0;tf): 

%  the  time  vector 
len=length(t); 
l=ones(1  Jen); 

%  Initialize  storage  for  Kalman  gains  and  the  states 
g=zeros(3Jen): 

X=zeros(3Jen); 

%  Emitter  location 

xt=input(’  Enter  the  x  coordinate  of  the  phone  in  km..’); 
yt=input(’  Enter  the  y  coordinate  of  the  phone  in  km..’); 
d=input(’Enter  the  delay  of  the  phone  in  microseconds..’); 
d=d*1e-6; 

xt=xt*1000; 

yt=yt*1000; 

%  A  priori  estimate  of  the  location  of  the  target  and  delay 
X0=[8000;8000;5e-6]; 

%  Location  of  the  cell  site  antennae 

xc=1 0000*1; 
yc=1 0000*1; 

%  Location  of  the  van  and  UAV 

xvan=10000*l+0*t; 

yvan=0*l+0*t; 
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dtheta=pi*rand(1  ,len)/1 .5; 
%dtheta=pi*.1; 

theta(1)=-pi/2; 

for  i=2:len 

theta(i)=theta(i-1  )+dtheta(i): 
end 

Cx=ioooo: 

Cy=10000; 

r=8000; 


xuav=Cx*l+(r*l).*cos(theta); 

yuav=Cy*l+(r*l).*sin(theta); 

%xuav=1 5*1+1  *t; 

%yuav=0*l; 

%  Calculate  the  TDOA  observations  for  the  two  receivers 

Z1  =(1/c)*(sqrt((xt-xc).^2+(yt-yc).^2)... 

+sqrt((xt-xvan).''2+(yt-yvan).‘^2)... 

-sqrt((xvan-xc).'^2+(yvan-yc).'^2))+d; 

Z2=(1/c)*(sqrt((xt*xc).-^2+(yt-yc).'^2)... 

+sqrt((xt-xuav).'^2+(yt-yuav).'^2)... 

-sq  rt((xuav-xc).^2+(yuav-yc)  .^2))+d; 

%  Inject  zero  mean  noise  into  the  observations 

n1  =sqrt(R1)*randn(1  ,len); 
n2=sqrt(R1  )*randn(1  Jen); 


Z1n=Z1+n1; 

Z2n=Z2+n2; 

%  Initialize  the  Kalman  Filter 

X(: ,  1  )=X0;  %  A  priori  estimate  of  the  emitter  location 

PK=PO; 

%  Initial  error  covariance  matrix 
%  calculate  the  filtered  estimates  of  the  emitter  location 
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for  l<=1  :len-1 , 

%  perform  the  first  calculation  for  the  measurement  at  the  van. 

Rtc=sqrt((X(1  ,k)-xc(k)).A2+(X(2,k)-yc(k)).'^2): 

Rtr=sqrt((X(1  ,k)-xvan(k)).^2+(X(2,k)-yvan(k)).'^2): 
Rrc=sqrt((xvan(k)-xc(k)).^2+(yvan(k)-yc(k)).'^2); 

h1x=(1/c)*((X(1,k)-xc(k))./Rtc+(X(1,k)-xvan(k))./Rtr-... 

(xvan(k)-xc(k))./Rrc); 

h1y=(1/c)*((X(2,k)-yc(k))./Rtc-(X(2,k)-yvan(k))./Rtr-... 

(yvan(k)-yc(k))./Rrc): 

HK=[h1xh1y  1]; 

%  calculate  the  estimate  of  TDOA  based  upon  the  estimate  of 
%  phone  location  and  location  of  van  and  UAV 

Zhat1  =(1  /c)*(sqrt((X(1  ,k)-xc(k)).^2+(X(2,k)-yc(k)).^2)... 
+sqrt((X(1,k)-xvan(k)).^2+(X(2,k)-yvan(k)).^2)... 
-sqrt((xvan(k)-xc(k)).^2+(yvan(k)-yc(k)).^2))+X(3,k); 

%  The  new  error  covariance  is  the  same  as  the  old  plus  Q 
PK=PK+Q; 

%  Calculate  the  Kalman  gains 
GK=PK*HK’*inv(HK*PK*HK’+R1 ); 

%  Calculate  the  updated  error  covariance  matrix 
%PK=(eye(2)-GK*HK)*PK: 

%  Joseph  form  of  covariance  update. 
PK=(eye(3)-GK*HK)*PK*(eye(3)-GK*HK)’+GK*R1*GK’; 

%  Calculate  the  smoothed  estimate  of  the  bearing  to  the  phone 
X(;,k)=X(:,k)+GK*(Z1  n(k)-Zhat1 ); 

X(:,k)=X(:.k); 

%  perform  the  calculations  for  the  measurement  at  the  UAV 

Rtc=sqrt((X(1,k)-xc(k)).A2+(X(2,k)-yc(k)).A2); 

Rtr=sqrt((X(1,k)-xuav(k)).'^2+(X(2,k)-yuav(k)).'^2); 

Rrc=sqrt((xuav(k)-xc(k)).'^2+(yuav(k)-yc(k)).'^2); 

h2x=(1/c)*((X(1  ,k)-xc(k))./Rtc+(X(1  ,k)-xuav(k))./Rtr-... 
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(xuav(k)-xc(k)),/Rrc); 

h2y=(1/c)*((X(2,k)-yc(k))./Rtc-(X(2,k)-yuav(k))./Rtr-... 

(yuav(k)-yc(k))./Rrc); 

HK=[h2x  h2y  1]; 

%  calculate  the  estimate  of  TDOA  based  upon  the  estimate  of 
%  phone  location  and  location  of  van  and  UAV 

Zhat2=(1/c)*(sqrt((X(1,k)-xc(k)).^2+{X(2,k)-yc(k)).^2)... 

+sqrt((X(1,k)-xuav(k)).''2+(X(2,k)-yuav(k)).'^2)... 

-sqrt((xuav(k)-xc(k))A2+(yuav(k)-yc(k))A2))+X(3,k); 

%  The  new  error  covariance  is  the  same  as  the  old  plus  Q 
PK=PK+Q; 

%  Calculate  the  Kalman  gains 
GK=PK*HK’*inv(HK*PK*HK’+R1 ); 

%  Calculate  the  updated  error  covariance  matrix 
%PK=(eye(2)-GK*HK)*PK: 

%  Joseph  form  of  covariance  update. 
PK=(eye(3)-GK*HK)*PK*(eye(3)-GK*HK)’+GK*RrGK’: 

%  Calculate  the  smoothed  estimate  of  the  bearing  to  the  phone 
X(:,k)=X(:,k)+GK*(Z2n(k)-Zhat2); 

X(:,k+1)=X(:,k); 

if  k==1 

Zhat=[Zhat1  ;Zhat2]; 

0lS6 

Zhat=[Zhat;Zhat1  ;Zhat2]; 
end 

%  save  the  error  covariance  matrices 
if  k=1 

P=PK(1:2,1:2); 

6lS6 

P=[P:PK(1:2,1:2)]; 

end; 

%  save  the  Kalman  gains 
g(:.k)=GK; 


no 


end; 


%  Plot  the  output 

Figure(1) 

cig: 

subplot(2,1,1) 

plot(t,X(1  ,:)/1 000,t,X(2,:)/1 000,t,xt*l/1 000,t.yt*l/1 000) 
title(’ Estimate  of  the  x  and  y  coordinates’): 
xlabel(’Number  of  observations’) 
ylabel(’Phone  location  in  km’) 

subplot(2,1,2); 
plot(t,X(3,:)/1e-6,t,d*l/1e-6) 
title(’ Estimate  of  the  phone  delay’); 
xlabel(’Number  of  obsen/ations’) 
ylabel(’delay  in  microseconds’) 
pause 

%  Steady  state  estimation  of  emitter  location 
xtavg=mean(X(1  ,ren-10:len)); 
ytavg=mean(X(2,len-1 0:len)); 

%  Plot  the  grid 
Figure(2) 

plot(xt/1000,yt/1000,’*’) 
hold  on 

plot(xvan(k)/1 000,yvan(k)/1 000, ’r*’) 
plot(X0(1  )/1 000,X0(2)/1 000,’+’) 


for  k=1:len-1; 

plot(xuav(k)/1 000,yuav(k)/1 000,’g*’); 

P1=[P(2*k-1,1)  P(2*k-1,2);P(2*k.1)  P(2*k,2)]: 
[xg,yg,Emax]=ellipa(P1 ,3,X(1  ,k),X(2,k));. 
plot(xg/1000,yg/1000,’b-’); 

if  rerri(k,15)==0 

[evan]=ellip1  (xc(k),yc(k),xvan(k),yvan(k),Zhat(2*k-1  ),c); 
[euav]=ellip1(xc(k),yc(k),xuav(k),yuav(k),Zhat(2*k),c); 
plot(evan(1  ,:)/1 000,evan(2,:)/1 000, ’g:’); 
plot(euav(1  ,:)/1 000,euav(2,:)/1 000,’y:’); 
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end 

end 

plot(X(1,:)./1000,X(2,:)./1000,’r’): 

axis(’equar) 

hold  off 

title(’Estimate  of  phone’) 

xlabel(’(km)’):ylabel(’(km)’); 

pause 

%  Plot  close  up  of  the  estimate  of  the  emitter  location 

Figure(3) 

cig 

axis([(xt- 1000)71000  (xt+1 000)71 000  (yt-1 000)71 000  ... 

(yt:i-1 000)71 000]); 
hold  on 

plot(xt71000,yt71000,’*’) 
plot(X(1  ,:).71 000,X(2,:).71 000); 


k=len-2; 

p=P(2*k-1:2*k.:); 

[xg,yg,Emax]=ellipa(p,C,X(1  ,len-2),X(2,len-2)); 
plot(xg71 000,yg71 000, ’b-’); 

[xg,yg,Emax]=ellipa(P1  ,C,X(1  ,len-1  ),X(2,len-1 )); 
plot(xg71000,yg71000,’b-’); 
plot(xtavg71 000,ytavg71 000,’r+’); 
hold  off 

title(’Close  up  of  position  of  the  Phone’) 
xlabel(’(km)’);ylabel(’(km)’); 

TDOA3  -  THE  COMBINED  SAT  AND  BURST  FILTER 
%  Michael  P.  Fallon 

%  This  program  combines  the  two  previous  methods.  The  SAT  and  burst 
%  TDOA  methods  are  combined  into  one. 

clear 

c=3.00e+08;  %  speed  of  light 

R1  =1 .694e-1 5;  %  The  covariance  of  measurement  error 

Q=5e+02*eye(2);  %  Variance  of  Plant  excitation  noise 

P0=5e+10*eye(2);  %  Initial  error  covariance 
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%  Choose  to  calculate  the  error  covariance  of  3  sigma 
C=1.0: 

tf=50: 

t=(0:tf): 

%  the  time  vector 

len=length(t): 

l=ones(1,len); 

%  Initialize  storage  for  Kalman  gains  and  the  states 
g=zeros(2Jen): 

X=zeros(2,Ien); 

%  Emitter  location 

xt=input(’  Enter  the  x  coordinate  of  the  phone  in  km..’); 
yt=input(’  Enter  the  y  coordinate  of  the  phone  in  km..’); 
d=0; 

d=d*1e-6; 

xt=xt*1000; 

yt=yt*1000; 

%  A  priori  estimate  of  the  location  of  the  target 
X0=[9000;9000]; 

%  Location  of  the  cell  site  antennae 

xc=1 0000*1; 
yc=1 0000*1; 

%  Location  of  the  van  and  UAV 
xvan=1 0000*l+0*t; 
yvan=0*l+0*t; 

dtheta=pi*rand(1  ,len)/2; 
theta(1)=-pi/2; 

for  i=2:len 

theta(i)=theta(i-1  )+dtheta(i); 
end 

Cx=10000; 

Cy=10000; 
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r=10000; 


xuav=Cx*l+(r*l).*cos(theta): 

yuav=Cy*l+(r*l).*sin(theta); 

%xuav=1 5*1+1  *t; 

%yuav=0*l; 

%  Calculate  the  TDOA  observations  for  the  two  receivers.  Z  is  the 
%  TDOA  for  the  burst  mode  and  Z1  and  Z2  represent  TDOA  for  the  SAT 
%  TDOA  for  the  van  and  UAV  respectively. 

Z 1  =(  1  /c)*(sqrt((xt-xc)  .^2+(yt-yc)  .^2) . . . 

+sqrt((xt-xvan)  .'^2+(yt-yvan).^2).. . 

-sqrt((xvan-xc)  .^2+(yvan-yc)  .^2))+d ; 

Z2=(1  /c)*(sqrt((xt-xc)  .''2+(yt-yc).'^2) . . . 

+sqrt((xt-xuav).^2+(yt-yuav).^2)... 

-sqrt((xuav-xc)  .'^2+(y  uav-yc)  .'^2))+d ; 

Z=(1/c)*(sqrt((xt-xvan).^2+(yt-yvan).'^2)-... 

sqrt((xt-xuav).'^2+(yt-yuav).'^2)); 

%  Inject  zero  mean  noise  into  the  observations 

n1=sqrt(R1)*randn(1  ,len); 
n2=sqrt(R1)*randn(1  ,len); 
n=sqrt(R1  )*randn(1  ,len); 

Z1n=Z1+n1; 

Z2n=Z2+n2; 

Zn=Z+n; 

%  Initialize  the  Kalman  Filter 

X(:,1  )=:X0;  %  A  priori  estimate  of  the  emitter  location 

PK=P0; 

%  Initial  error  covariance  matrix 

%  calculate  the  filtered  estimates  of  the  emitter  location 
for  k=1  :len-1 , 

%  perform  the  first  calculation  for  the  measurement  at  the  van. 

Rtc=sqrt((X(1.k)-xc(k)).A2+(X(2,k)-yc(k)).A2); 

Rtr=sqrt((X(1,k)-xvan(k)).'^2+(X(2,k)-yvan(k)).^2): 

Rrc=sqrt((xvan(k)-xc(k)).''2+(yvan(k)-yc(k)).'^2); 


114 


h1  x=(1/c)*((X{1  ,k)-xc(k))./Rtc+(X(1  ,k)-xvan(k))./Rtr-... 
(xvan(k)-xc(k))./Rrc); 

h1y=(1/c)*((X(2,k)-yc(k))./Rtc-(X(2,k)-yvan(k))./Rtr-... 

(yvan(k)-yc(k))./Rrc): 

HK=[h1xh1y]; 

%  calculate  the  estimate  of  TDOA  based  upon  the  estimate  of 
%  phone  location  and  location  of  van  and  UAV 

Zhat1=(1/c)*(sqrt((X(1,k)-xc(k))A2+(X(2,k)-yc(k)).^2)... 

+sq  rt((X(  1 ,  k)-xvan  (k))  .'^2+(X(2,  k)-yvan  (k))  .-^2) . . . 
-sqrt((xvan(k)-xc(k)).'^2+(yvan(k)-yc(k)).^2)); 

%  The  new  error  covariance  is  the  same  as  the  old  plus  Q 
PK=PK+Q: 

%  Calculate  the  Kalman  gains 
GK=PK*HK’*inv(HK*PK*HK’+R1 ); 

%  Calculate  the  updated  error  covariance  matrix 
%PK=(eye(2)-GK*HK)*PK; 

%  Joseph  form  of  covariance  update. 
PK=(eye(2)-GK*HK)*PK*(eye(2)-GK*HK)’+GK*RrGK’; 

%  Calculate  the  smoothed  estimate  of  the  bearing  to  the  phone 
X(:,k)=X(:,k)+GK*(Z1  n(k)-Zhat1 ); 

X(:,k)=X(:,k); 

%  perform  the  calculations  for  the  measurement  at  the  UAV 

Rtc=sqrt((X(1  ,k)-xc(k)).A2+(X(2,k)-yc(k)).A2); 

Rtr=sqrt((X(1,k)-xuav(k)).^2+(X(2,k)-yuav(k)).'^2); 

Rrc=sqrt((xuav(k)-xc(k)).^2+(yuav(k)-yc(k)).'^2); 

h2x=(1/c)*((X(1  ,k)-xc(k))./Rtc+(X(1  ,k)-xuay(k))./Rtr-... 
(xuav(k)-xc(k))./Rrc); 

h2y=(1/c)*((X(2,k)-yc(k))./Rtc-(X(2,k)-yuav(k))./Rtr-... 

(yuav(k)-yc(k))./Rrc); 

HK=[h2x  h2y]; 

%  calculate  the  estimate  of  TDOA  based  upon  the  estimate  of 
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%  phone  location  and  location  of  van  and  UAV 


Zhat2=(1/c)*(sqrt((X(1,k)-xc(k)).^2+(X(2,k)-yc(k))A2)... 

+sqrt((X(  1 ,  k)-xuav(k))  A2+{X{2,  k)-yuav(k))  .^2) . . . 
-sqrt((xuav(k)-xc(k)).^2+(yuav(k)-yc(k)).^2)); 

%  The  new  error  covariance  is  the  same  as  the  old  plus  Q 
PK=PK+Q; 

%  Calculate  the  Kalman  gains 
GK=PK*HK’*inv(HK*PK*HK’+R1 ); 

%  Calculate  the  updated  error  covariance  matrix 
%PK=(eye(2)-GK*HK)*PK; 

%  Joseph  form  of  covariance  update. 
PK=(eye(2)-GK*HK)*PK*(eye(2)-GK*HK)’+GK*R1*GK’; 

%  Calculate  the  smoothed  estimate  of  the  bearing  to  the  phone 
X(:,k)=X(:,k)+GK*(Z2n(k)-Zhat2); 

X(:,k)=X(:,k); 

Rvan=sqrt((X(1,k)-xvan(k)).^2+(X(2,k)-yvan(k)).'^2): 

Ruav=sqrt((X(1,k)-xuav(k)).^2+(X(2,k)-yuav(k)).^2); 

hx=(1/c)*((X(1  ,k)-xvan(k))./Rvan-(X(1  ,k)-xuav(k))./Ruav); 
hy=(1/c)*((X(2,k)-yvan(k))./Rvan-(X(2,k)-yuav(k))./Ruav); 

HK=[hx  hy]; 

%  calculate  the  estimate  of  TDOA  based  upon  the  estimate  of 
%  phone  location  and  location  of  van  and  UAV 

Zhat=(1 /c)*(sqrt((X(1  ,k)-xvan(k))^2+(X(2,k)-yvan(k))^2)-. . . 

sqrt((X(1,k)-xuav(k))^2+(X(2,k)-yuav(k))^2)); 

Zh(k)=Zhat; 

%  The  new  error  covariance  is  the  same  as  the  old  plus  Q 
PK=PK+Q; 

%  Calculate  the  Kalman  gains 
GK=PK*HK’*inv(HK*PK*HK’+R1 ); 

%  Calculate  the  updated  error  covariance  matrix 
%PK=(eye(2)-GK*HK)*PK; 
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%  Joseph  form  of  covariance  update. 
PK=(eye(2)-GK*HK)*PK*(eye(2)-GK*HK)’+GK*RrGK’: 

%  Calculate  the  smoothed  estimate  of  the  bearing  to  the  phone 
X(:,k)=X(:,k)+GK*(Zn(k)-Zhat); 

X(:,k+1)=X(:,k); 


if  k==1 

Zhat=[Zhat1:Zhat2]; 

else 

Zhat=[Zhat;Zhat1  ;Zhat2]: 
end 

%  save  the  error  covariance  matrices 
if  k==1 
P=PK; 
else 

P=[P;PK]; 

end: 

%  save  the  Kalman  gains 
g(:.k)=GK; 


end; 

%  Plot  the  output 

Figure(1) 

cig; 

plot(t,X(1  ,:)/1 000,t,X(2,:)/1 000) 

title(’ Estimate  of  the  x  and  y  coordinates’); 

xlabel(’Number  of  observations’) 

ylabel(’Phone  location  in  km’) 

grid 

pause 

%  Steady  state  estimation  of  emitter  location 
xtavg=mean(X(1  ,len-10:len)); 
ytavg=mean(X(2,len-1 0:len)); 

%  Plot  the  grid 
Figure(2) 

plot(xt/1000,yt/1000,’*’) 
hold  on 

plot(xvan(k)/1 000,yvan(k)/1 000,’r*’) 
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plot(X0(1  )/1 000,X0(2)/1 000,’+’) 
for  k=1  :len-1 ; 

plot(xuav(k)/1 000,yuav(k)/1 000, ’g*’); 

P1=[P(2*k-1,1)  P(2*k-1,2);P(2*k,1)  P(2*k,2)]; 
[xg,yg,Emax]=eflipa(P1  ,C,X(1  ,k),X(2,k)); 
plot(xg/1000,yg/1000,’b-’); 

end 

plot(X(1  ,:)./1 000,X(2,:)./1 000,’r’); 

axis(’equar) 

hold  off 

title(’Estimate  of  phone’) 

xlabel(’(km)’);ylabel(’(km)’); 

pause 

%  Plot  close  up  of  the  estimate  of  the  emitter  location 

Figure(3) 

cig 

axis([xt/1000-.1  xt/1000+.1  yt/1000-.1  yt/1000+.1]) 
hold  on 

plot(xt/1000,yt/1000,’*’) 
plot(X(1  ,:)./1000,X(2,;)./1000); 


k=len-2; 

p=P(2*k-1:2*k,:); 

[xg,yg,Emax]=ellipa(p,C,X(1  ,len-2),X(2,len-2)); 
plot(xg/1 000,yg/1 000,’b-’); 

[xg,yg,Emax]=ellipa(P1  ,C,X(1  ,len-1),X(2,len-1)); 
plot(xg/1 000, yg/1 000,’b-’); 
plot(xtavg/1 000,ytavg/1 000,’r+’); 
hold  off 

title(’Close  up  of  position  of  the  Phone’) 
xlabel(’(km)’);ylabel(’(km)’); 

ELIPA  -  THE  ERROR  ELLIPSE  SUBROUTINE 

%  ellipa.m 
%  Michael  P.  Fallon 

%  this  program  calculates  the  error  ellipsoids  given 
%  the  error  covariance  matrix  ahd  the  location  of  the  estimate. 

%  This  method  is  based  upon  Kirk’s  work  using  eigenvalues  and 
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%  eigenvectors. 

% 

function[xout,yout,Emax]=ellipa(PK,c,xt,yt) 

[V,lam]=eig(PK); 
theta=atan(V(2,1  )/V(1 ,1 )); 
t=0:.01:2*pi; 

x=sqrt(lann(1 ,1  )*c^2)*cos(t-theta); 
y=sqrt(lam(2,2)*c^2)*sin(t-theta); 

xout=x+xt; 

yout=y+yt: 

Emax=max(abs(V(:,1))); 

ellipl  -  the  position  ellipse  subroutine 

%  ellipl  .m 
% 

%  This  function  plots  the  possible  positions  for  the  phone  based  on 
%  the  geometry  relating  to  the  SAT  channel  method. 

% 

%  cx.xy  position  of  the  cell  antennae 

%  rx,ry  position  of  the  receiver  (van  or  UAV) 

%  zhat  TDOA 

%c  speed  of  light 
%  X  The  position  of  the  ellipse. 

% 

function[x]=ellip1(cx,cy,rx,ry,zhat,c) 

cenb(=(cx+rx)/2; 

centy=(cy+ry)/2; 

d=sqrt((cx-rx)^2+(cy-ry)^2); 

a=(c*zhat+d)/2; 

e=d/2; 

b=sqrt(a'^2-e^2); 

phi=0;pi/100:2*pi; 

n=sin(phi); 

n=n.'^2; 

n=n*a''2; 

m=(b^2)*((cos(phi)).'^2): 
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r=(aA2*bA2)./(n+m): 

r=sqrt(r); 

%r=sqrt(a'^2*b^2/(a'^2.*(sin(phi)).^2+b^2.*(cos(phi)).^2)): 


x=r.*cos(phi); 

y=r.*sin(phi): 

x1=x+centx; 
y1  =y+centy; 

x=[x1:y1]; 
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